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1.  INTRODUCTION 


The  need  to  develop  a  technology  base  for  producing  damage  tolerant  composite  materials 
is  viewed  as  a  key  issue  in  obtaining  wide  acceptance  of  advanced  composites  for  aircraft  and 
aerospace  structural  applications.  This  need  is  further  emphasized  by  the  damage  tolerant 
requirements  set  forth  which  must  be  met  with  if  the  material  is  to  be  used  in  primary  structures. 
Research  work  in  recent  years  has  demonstrated  that  failure  processes  associated  with  delamination 
are  of  primary  concern  in  structural  composite  materials.  In  order  to  predict  the  onset  of 
delamination  and/or  propagation  of  delamination  in  complex  geometries  (e.g.  around  a  crack  or 
circular  hole),  a  three-dimensional  stress  analysis  model  is  required  which  can  be  applied  to  thick 
laminates  (i.e.,  laminates  containing  more  than  six  plies).  Although  finite  element  codes  based  on 
three-  dimensional  codes  are  available,  they  are  very  complex,  time  consuming  and  costly.  In 
addition,  many  of  these  codes  are  not  capable  of  analyzing  thick  laminates.  The  global-local  model 
described  in  this  report  has  that  capability.  To  date,  however,  the  solutions  have  been  limited  to 
simple  straight  edges  of  laminates.  In  order  to  fully  utilize  this  approach,  numerical  procedures  are 
required  to  solve  problems  with  curved  boundaries  and  sharp  discontinuities. 

Another  key  issue  in  today's  composite  technology,  particularly  in  the  utilization  of 
carbon-carbon  composites  in  exit  cone  applications,  is  the  p  rocessing.  There  is  a  history  of  failures 
in  three-dimensional  carbon-carbon  exit  cones  during  processing.  Such  failures  are  not  predictable 
by  current  models.  Two  classes  of  models  are  required  in  the  analytical  treatment  of  composite 
bodies  under  processing  conditions.  The  first  involves  the  study  of  transport  phenomena.  The 
second  is  a  mechanical  model  which  can  examine  the  stress,  strain  and  displacement  fields  within  the 
body  and  forms  the  basis  for  failure  prediction.  Both  these  models  are  critical  in  determining  the 
effect  of  processing  parameters  on  resulting  residual  stresses  and  developing  failure  criteria. 

In  this  report,  a  description  of  research  efforts  towards  the  objective  of  developing 
three-dimensional  models  capable  of  describing  the  behavior  of  thick  laminated  composites  as  well  as 
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the  residual  stresses  in  carbon-carbon  composites,  is  given.  On  a  broader  level,  both  macro-  and 
micro-mechanics  approaches  have  been  considered.  The  macro-mechanics  approach  deals  with  the 
problem  of  stress  analysis  in  composite  laminates  under  different  loading  conditions  with  various 
types  of  boundary  conditions  (including  complex  ones  such  as  stress-free  edges).  The  solutions  to 
these  problems  are  attempted  by  using  numerical  techniques  such  as  finite  difference  and  finite 
element  methods.  The  micro-mechanics  approach  deals  with  the  determination  of  effective 
thermoelastic  properties  and  the  stress  distribution  under  three-dimensional  mechanical  and/or 
hygrothermal  loading. 

The  macro-mechanics  study  in  this  report  includes  an  analytical  procedure  leading  to  a 
qualitative  understanding  of  the  nature  of  stress  field  near  the  regions  of  steep  stress  gradients.  One 
aspect  of  this  study  is  the  free  edge  problem  in  a  multi-layered  composite  laminate.  The  potentially 
high  stress  gradients  near  the  free  edges  may  limit  the  load  carrying  capability  of  the  structure.  It 
may  also  become  a  source  of  laminate  failure.  The  presence  of  complex  stress  fields  near  the  free 
edges  of  laminates  has  been  established  experimentally.  However,  the  analytical  investigation  of 
the  nature  and  magnitude  of  the  free  edge  stresses  by  various  researchers  has  lead  to  disparate 
conclusions.  As  an  extension,  a  reliable  analytical  procedure  for  solving  the  free  edge  problem  may 
also  lead  to  an  understanding  of  the  various  failure  modes  that  have  been  shown  to  occur  in 
composite  laminates  [1,2,3,4].  Also,  in  practical  applications  of  composite  materials,  a  large 
number  of  anisotropic  layers  are  usually  present.  Consequently,  an  exact,  three-dimensional 
analysis  of  such  a  laminate  becomes  extremely  complicated.  The  references  mentioned  before 
indicate  the  importance  of  defining  the  stress  fields  in  various  layers  of  multi-layered  laminate. 
Hence  it  becomes  important  to  establish  an  analytical  procedure  to  consider  all  these  aspects  of 
laminate  behavior.  In  order  to  obtain  an  understanding  regarding  the  p’T)cedures  that  are  available 
for  solving  free  edge  problems  of  laminates,  a  brief  review  or  literature  prertaining  to  general  laminate 
stress  analysis  is  given. 

The  basis  for  the  micro-mechanics  study  in  this  report  is  the  fact  that  the  mechanical 
properties  of  a  fiber-reinforced  material,  including  the  failure  modes  and  mechanisms,  are  governed 
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in  part  by  the  transfer  of  stresss  between  fiber  and  the  matrix.  This  transfer  occurs  across  the 
interface  between  the  components,  and  the  properties  of  this  interface,  therefore,  will  affect  the 
properties  of  the  composite.  For  example  the  strength  of  the  interface  in  tension  has  a  direct  bearing 
on  the  composite  transverse  and  compressive  strengths  and  delamination  parallel  to  the  fibers; 
whereas,  the  interface  shear  strength  principally  affects  the  load  transfer  length;  composite  fracture 
under  conditions  of  fiber  pullout  and  the  deformation  of  the  matrix.  Recently,  there  has  been  an 
increasing  use  of  coated  fibers  as  a  reinforcement  in  some  new  application  areas  such  as 
carbon-carbon  composites,  electric  composites,  metal-matrix  composites  or  ceramic  composites 
intended  for  high  temperature  applications.  The  coating  applied  on  the  fiber  serves  to  enhance  the 
electrical  conductivity,  or  the  elastic  stiffness  or  strength,  or  simply  to  serve  as  a  chemical  reaction 
barrier.  From  this  point  of  view,  it  is  essential  to  understand  the  interaction  between  the  coating, 
fiber  and  matrix  in  order  to  arrive  at  conlusions  regarding  the  nature  of  stress  distribution  and  the 
consequent  failure  mechanisms  under  mechanical  and  hygrothermal  laods.  This  report  includes  a 
model  to  study  this  aspect  of  composite  behavior. 

1.1  LITERATURE  REVIEW 


There  exist  in  literature  approximate  theories  for  achieving  rea.sonably  accurate  laminate 
stress  analysis.  The  most  popular  of  these  is  the  classical  laminate  theory  which  is  an  extension  of 
the  well-known  classical  plate  theory.  This  procedure  [5],  which  included  the  bending-stretching 
coupling,  later  incorporated  laminate  shear  deformation  thus  resulting  in  a  generalized 
Reissner-Mindlin  plate  theory  [6].  In  these  theories,  the  laminate  was  assumed  to  be  in  a  state  of 
plane  stress.  The  assumptions  involved  in  these  theories  were  found  to  be  too  restrictive  for  general 
application  [7,8].  They  lead  to  inaccurate  prediction  of  interlaminar  stresses  at  the  free  edges, 
although  they  yield  reasonbly  accurate  results  for  regions  away  from  the  free  edges.  A  higher  order 
plate  theory  [9]  in  which  the  displacements  no  longer  vary  linearly,  was  applied  [10]  to  evaluate  the 
interlaminar  normal  stress  distribution  in  a  free  edge  boundary  value  problem,  but  only  on  a  plane  of 
symmetry.  In  the  above  mentioned  theories,  a  displacement  field  is  generally  assumed  that  is 
continuous  across  the  entire  thickness  of  the  laminate.  This,  however,  does  not  guarantee  continuity 
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of  tractions  at  the  interfaces  between  layers.  Also,  the  traction  boundary  conditions  for  these 
approaches  are,  in  general,  shown  to  be  insufficient  to  guarantee  the  equilibrium  of  edge  boundary 
subregions  [10]. 

The  study  of  delamination  phenomenon  in  structural  composite  laminates  began  with 
analytical  and  experimental  observations  of  the  response  of  such  bodies  in  the  vicinity  of  a  free 
edge.  In  1967,  Hayashi  [11]  presented  the  first  analytical  model  treating  interlaminar  stresses  in 
what  has  come  to  be  known  as  the  "free  edge  problem".  Characteristically,  this  work  focused  on  the 
computation  of  interlaminar  shear  stress,  as  in  the  early  stages  of  composite  research  interlaminar  and 
delamination  effects  were  viewed  as  being  synonymous  with  interlaminar  shear.  The  presence  of 
interlaminar  normal  stress,  being  of  a  more  subtle  origin  and  also  seemingly  defying  common 
intuition,  was  not  appreciated  until  many  years  after  the  pioneering  work  of  Hayashi.  The 
development  of  the  Hayashi  model  was  based  upon  the  implicit  assumption  that  the  in-plane  stresses 
within  a  given  layer  did  not  depend  upon  the  thickness  coordinate.  The  magnitude  of  the  maximum 
interlaminar  shear  stress  was  calculated  to  be  a  relatively  large  value  in  a  glass  epoxy  [0/90]s 
laminate.  Unfortunately,  however,  owing  to  the  omission  of  the  interlaminar  normal  stress,  the 
computed  stress  field  within  each  layer  does  not  satisfy  moment  equilibrium. 

The  first  reported  experimental  observations  involving  free  edge  delamination  were  made 
by  Foye  and  Baker  [12].  In  that  work,  tremendous  differences  in  fatigue  life  of  boron-epoxy 
composite  laminates  as  a  function  of  layer  stacking  sequence  were  reported.  Severe  delaminations 
were  witnessed  in  that  work  and  were  identified  as  the  primary  source  of  strength  degradation  in 
fatigue. 


With  regard  to  free  edge  problems,  a  finite  width  laminate  under  uniform  extensional 
strain  was  first  analyzed  by  Pipes  and  Pagano  [13].  A  finite  difference  scheme  was  used  based  on  a 
formulation  assuming  isotropic  elasticity.  Subsequently,  various  authors  have  conducted 
investigations  on  this  problem  by  using  various  techniques  such  as  finite  difference,  finite  element 
and  series  solutions.  References  [2,  3, 4, 14]  lead  to  an  understanding  of  various  failure  modes  that 
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have  been  shown  to  occur  in  composite  laminates.  The  bending  and  stretching  of  laminate 
deformation  and  Reissner-Mindlin  plate  theory  have  been  incorporated  [5,  6]  and  higher  order 
classical  laminate  theory  [9, 10]  was  applied  to  evaluate  the  interlaminar  normal  stress  distribution, 
but  only  at  a  plane  of  symmetry.  When  a  large  number  of  layers  is  present  in  the  laminate 
construction,  contemporary  models  are  incapable  of  providing  precise  solution  of  the  local  stress 
fields  in  the  vicinity  of  free  edge.  In  an  effective  moduli  global  model  [14],  only  the  extensional 
response  of  the  regions  was  conducted,  i.e.  the  flexural  and  flexural-extension  coupling 
characteristics  of  laminated  bodies  were  ignored.  A  global  representation  of  a  three  dimensional 
laminate  model  [15]  is  a  generali2ation  and  improvement  of  the  material  model  given  in  reference  14. 
In  the  literature  on  the  finite  element  solutions  [14, 16],  the  traction  free  edge  conditions  are  satisfied 
in  weighted  integrated  sense  or  lead  to  oscillating  solutions  near  the  edges.  Singular  hybrid  stress 
finite  element  models  are  also  employed  [17, 18, 19].  The  free  edge  stresses  in  layered  plates  have 
also  been  evaluated  by  using  eight  nodes  isoparametric  elements  [20],  It  was  shown  [19,  20]  that 
the  laminate  idealization  for  a  reasonably  accurate  finite  element  analysis  had  to  be  very  fine,  i.e.  a 
quarter  of  the  laminate  was  divided  into  about  600  elements.  No  more  than  six  layers  would  be 
considered  for  numerical  calculations.  For  moderately  large  number  of  plies,  these  approaches  will 
lead  to  computer  storage/economic  difficulties.  In  other  studies  a  singular  hybrid  finite  element 
model  [21]  is  employed.  The  finite  difference  method  reported  [13]  leads  to  difficulties  in  calculating 
transverse  shear  stresses  at  the  free  edges.  More  recently,  a  three-dimensional  finite  difference 
scheme  has  also  been  developed  [22].  In  addition,  a  perturbation  technique  [23]  and  a  series 
solution  [24]  have  also  been  presented  for  obtaining  the  solution  to  the  free  edge  problem. 

In  order  to  circumvent  the  complexity  of  an  exact  three  dimensional  elastic  analysis,  while 
at  the  same  time  being  a  reasonably  precise  method  of  studying  the  stress  fields  in  laminate  with 
moderately  large  number  of  layers,  a  unified  tractable  model  for  the  elastic  response  of  the  individual 
plies  of  the  laminate  has  been  introduced  [24].  But  this  model,  termed  the  local  model,  in  which 
each  layer  is  represented  as  a  homogeneous,  anisotropic  continuum,  becomes  intractable  as  the 
number  of  layers  becomes  large.  To  overcome  this  difficulty,  the  local  model  has  been  recently 
extended  into  a  formulation  of  global-local  variational  model  [25].  In  this  model,  a  predetermined 
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area,  termed  local  region,  of  interest  is  represented  by  Reissner's  functional  that  involves  both 
stresses  and  displacements.  In  this  local  model,  the  interface  continuity,  including  three 
displacement  components  and  the  nonnal  and  shear  stress  components  associated  with  the  thickness 
direction  are  satisfied.  For  a  large  number  of  layers  (exceeding  6),  the  exponential  parameter 
becomes  very  large  and  consequently  causes  computer  overflow. 

On  the  other  hand,  for  the  rest  of  the  region  other  than  the  local  region,  termed  global 
region,  the  principle  of  potential  energy  is  applied  following  an  assumed,  elementary  displacement 
field  and  a  consideration  of  effective  or  smeared  laminate  moduli.  While  the  local  model  employs 
the  theory  presented  in  [24],  the  global  model  is  based  on  higher  order  laminate  theory  given  in  [9]. 
The  present  report  addresses  the  problem  encountered  by  the  previous  methods  of  reasonably 
accurate  analysis  of  composite  laminates  [26].  In  references  25  and  26,  a  class  of  boundary  value 
problems  with  free  edges  has  been  solved  and  the  results  are  presented.  This  numerical  problem 
pertains  to  the  case  of  a  laminate  of  finite  width,  similar  to  that  of  the  case  of  plane  strain,  in  which 
the  stresses  depend  only  on  the  width  and  thickness  coordinates.  The  laminate  is  subjected  to  a 
uniform  axial  strain  in  the  length  direction.  The  solution  was  sought  in  the  form  of  an  exponential 
series  for  the  field  variables.  It  was  found  [25]  that  for  a  large  number  of  layers,  the  exponential 
parameter  becomes  very  large  with  a  consequent  exceeding  of  the  computer  limits.  In  other  words, 
the  number  of  layers  in  the  laminate  had  to  be  restricted.  The  largest  number  of  layers  that  could  be 
considered  by  this  numerical  solution  technique  was  6. 

This  report  describes  some  attempts  at  alternative  solution  techniques  to  circumvent  the 
above  mentioned  difficulty.  One  of  the  techniques  is  based  on  the  use  of  an  available  computer 
library  subroutine  to  solve  a  set  of  linear  differential  equations  with  an  appropriate  number  of 
boundary  conditions.  Yet  another  method  involves  the  finite  difference  modeling  of  the 
one-dimensional  problem  mentioned  before.  All  these  methods  are  applicable  to  the  local  model 
only  when  each  layer  of  the  laminate  is  considered  individually  with  regard  to  the  satisfaction  of  the 
layer  equilibrium  and  the  interlayer  continuity  conditions.  All  these  are  satisfied  exactly  once  the  set 
of  in-plane  stress  components  are  assumed  and  Reissner's  variational  functional  invoked. 
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The  present  report  also  describes  a  numerical  method  of  analysis  which  is  applicable  to 
the  global-local  model  of  the  laminate.  This  method  involves  a  discretization  process  using  finite 
element  technique.  A  current  research  literature  survey  was  conducted  with  particular  reference  to 
the  finite  element  models  based  on  mixed  variational  principles  (Reissner's  principle,  in  this  ca.se). 
The  purpose  of  this  literature  review  was  to  establish  the  guidelines  for  the  application  of  present 
finite  element  technique  to  the  global-local  model  of  the  laminate. 

It  has  been  established  that  the  failure  mechanisms  and  the  mechanical  properties  of 
fiber-reinforced  composites  are  controlled  in  part  by  the  transfer  of  stress  between  the  fiber  and  the 
matrix  [27].  For  theoretical  analysis  (a  micro-mechanics  approach),  the  interphase  region  between 
the  fiber  and  the  matrix  can  be  modeled  as  a  c(»ting  [28].  A  single  fiber  test  [29]  has  frequently  been 
used  to  characterize  the  fiber-matrix  interface  (  or  rather '  interphase  ’  since  the  region  adjacent  to 
the  fiber  has  its  unique  properties ).  Quite  often,  the  interphase  region  is  a  product  of  the  processing 
conditions  involved  in  the  manufacture  of  the  composite. 

In  the  case  of  continuous  fiber  composite,  a  number  of  works  [30  -34]  have  been  made  to 
compute  the  stress  field  in  a  composite  subjected  to  thermo-mechanical  loadings,  or  to  predict  its 
stiffness.  The  model  used  in  the  above  works  is  either  a  two-phase  model  consisting  of  inner 
cylinder  with  fiber  properties  and  an  outer  cylinder  with  the  properties  of  the  matrix  ,  or  a 
three-phase  model  in  which  one  more  cylinder  is  added  to  the  outside  of  the  two-phase  model  with 
the  composite  properties.  Recently,  Mikata  and  Taya  [35]  have  studied  the  stress  field  in  a  coated 
continuous  fiber  composite  which  requires  a  four-phase  model  that  consists  of  fiber,  coating,  matrix 
and  surrounding  composite  body.  The  material  properties  of  the  surrounding  body  (composite)  were 
obtained  by  using  a  rule  of  mixtures.  The  solution  to  the  stress  distribution  was  determined  with  the 
composite  subjected  to  three  independent  boundary  conditions,  namely,  axisymmetric  temperature 
change,  uniaxial  applied  stress  and  biaxial  applied  stress. 
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1.2  PROJECT  SUMMARY 


This  report  contains  analytical  models  developed  for  the  prediction  of  three-dimensional 
response  of  composites.  Two  aspects  of  the  problem  of  stress  and  failure  analysis  are  considered. 
One  pertains  to  a  ntacromechanical  model  in  which  the  development  of  a  global-local  finite  element 
method  for  stress  analysis  of  composite  laminates  is  given.  The  other  aspect  incorporates  a  micro¬ 
mechanics  approach  leading  to  a  model,  termed  NDSANDS,  for  the  stress  and  failure  analysis  of 
carbon-carbon  composites.  Again,  the  global-local  f(»mulation  for  the  laminate  stress  analysis  is 
numerically  solved  by  two  different  finite  element  techniques.  The  NDSANDS  model  incoiporates 
the  effects  of  prcx;essing  by  evaluating  the  residual  stresses  due  to  hygrothermal  effects.  A  number  of 
problems  has  been  solved  for  illustrating  the  effectiveness  all  the  above  mentioned  models.  It  is 
believed  that  these  models  are  of  extreme  importance  in  studying  the  optimum  design  parameters  in 
structural  composites. 
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2.  SOLUTION  TO  LOCAL  MODEL 


This  section  refers  to  the  attempts  for  obtaining  solution  to  the  local  model  [24]. 
Specifically,  the  first  part  of  this  section  deals  with  an  attempt  to  obtain  a  direct  solution  by  using  a 
library  subroutine  to  solve  a  set  of  linear  equations.  The  second  part  refers  to  a  one-dimensional 
finite  difference  scheme  to  solve  the  problem  of  an  infinitely  long  laminate  subjected  to  uniform  axial 
strain.  As  mentioned  in  the  previous  section,  the  exponential  series  solution  to  the  field  equations 
for  the  local  model  of  the  laminate  encountered  difficulties  when  implemented  on  the  computer.  This 
has  been  reported  in  reference  26  where  a  class  of  boundary  value  problems  with  free  edges  has  been 
solved  and  the  results  presented.  This  numerical  problem  pertains  to  the  case  of  a  laminate  of  finite 
width,  similar  to  that  of  the  case  of  plane  strain,  in  which  the  stresses  depend  only  on  the  width 
coordinate  y  and  the  thickness  coordinate  z.  The  laminate  is  subjected  to  a  uniform  axial  strain  in  the 
x-direction.  The  solution  was  sought  in  the  form  of  an  exponential  series  as 

fk  =  F*  e^y  (2.1) 

where  fk  represents  the  dependent  variables  for  the  k  th  layer  and  pk  denotes  the  coefficients  for  the 
corresponding  layer.  After  substituting  eq.  (2.1)  into  the  appropriate  field  equations  obtained  from 
theoretical  developments  for  the  local  model  [24]  the  values  of  X  were  obtained  by  setting  the 
determinant  of  the  coefficients  to  zero.  It  was  found  [24]  that  for  a  large  number  of  layers,  X  became 
so  very  large  as  to  exceed  the  computer  limits.  In  other  words,  the  number  of  layers  in  the  laminate 
had  to  be  restricted.  The  largest  number  of  layCTS  that  could  be  considered  by  this  numerical  solution 
technique  was  six.  The  present  section  discusses  an  alternative  mathematical  approach  to  the 
solution  to  circumvent  the  above  mentioned  difficulty.  This  method  invloves  the  use  of  an  available 
computer  library  subroutine  to  solve  a  set  of  linear  differential  equations. 
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2.1  DIRECT  SOLUTION 


One  of  the  alternative  solution  procedures  (direct  solution)  considered  makes  use  of  an 
IMSL  subroutine,  DVCPR,  to  solve  a  s^  of  equadtms.  The  purpose  of  this  subroutine  is  to  solve  a 
system  of  ordinary  differential  equations  with  boundary  conditions  defined  at  two  points.  The 
particular  sulnoutine  DVCPR,  for  example,  is  based  on  a  variable  step  size  finite  difference  method. 
The  fundamental  idea  was  to  cast  the  set  of  second  order  linear  differential  equations  governing  the 
behavior  of  each  layer  of  the  laminate  into  a  set  of  single  order  equations.  This  process  would  result 
in  a  set  of  additional  equations  and  all  the  equations  could  then  be  solved  simultaneously  by  using  the 
computer  subroutine. 

The  variables  in  the  set  of  governing  equations  are  the  weighted  displacement  functions 
and  the  interlaminar  stress  components .  For  any  layer,  these  are  U,  V,  W,  Q ,  y ,  X  •  <)>  >  Pi  •  P2> 
$1  <  S2 ,  ti ,  t2  where  the  Erst  six  refer  to  the  weighted  displacement  functions  and  the  three  pairs  of 
the  last  six  refer  to  the  interlaminar  stresses  associated  with  the  thickness  direction.  The  numerical 
problem  for  solution  consists  of  a  finite  width  laminate  subjected  to  a  uniform  axial  strain  =  e. 
By  virtue  of  the  geometry,  the  stress  components  and  the  displacements  depend  only  on  y  and  z 
coordinates.  Consequently,  all  the  derivatives  in  the  governing  equations  are  with  respect  to  y 
only,  thus  reducing  the  solution  procedure  to  that  of  a  one-dimensional  problem.  It  can  be  seen 
from  equations  (6)  through  (9)  of  reference  26  that  there  exist  single  order  derivatives  in  the 
interlayer  continuity  conditions  and  the  boundary  conditions  on  the  upper  and  lower  surfaces. 
However,  in  the  layer  equilibrium  equations  (5)  of  ref«ence  26,  the  highest  order  of  derivative  is  2. 
Hence,  a  new  set  of  variables  was  Erst  written  as,  for  example, 

U  =  U’ 

V  =  V’  etc.  (2.1.1) 

where  prime  refers  to  differentiation  with  respect  to  y.  The  new  variables  with  bars  on  them  were 
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also  considered  as  unknowns  in  further  formulation.  Thus,  the  layer  equilibrium  equations  were 
written  as  single  order  differential  equations  in  terms  of  variables  U,  U  etc. 


In  accordance  with  the  requirements  of  the  computer  subroutine,  the  next  step  was  to 
derive  a  set  of  equations  of  the  form: 


(Uk)'  =  f(Uk,Vk.Wk.Qk,  x|/k ,  Uk,^,Wk, 

Pl*^.P2^Sik,S2k,tik,t2k  )  (2.2) 


where  superscript  k  refers  to  the  layer. 

It  is  to  be  noted  that  the  right  hand  side  functicm  in  (2.2)  does  not  contain  any  derivative. 

Similar  equations  can  be  written  for  each  of  the  dependent  variable  appearing  with  function 
parantheses  in  the  above  equation.  It  was  found  during  the  course  of  such  a  procedure  that  one  of 
the  weighted  displacement  functions,  viz.,  <|>  could  be  eliminated  by  using  the  lower  surface 
condition  ( fw  obtaining  <{>1 )  and  a  judicious  combination  of  intersurface  continuity  conditions  and 
the  particular  layer  equilibrium  equation  that  involves  (for  obtaining  <|>k,  k  >  1 ).  Thus,  only 
six  of  the  seven  layer  equilibrium  equations  needed  to  be  considered  for  final  solution.  The  <}>k's 
obtained  in  this  manner  were  used  elsewhere. 

Also,  two  more  sets  of  equations  were  obtained,  one  by  equating  the  expression  for 
(<|>k)’  from  the  interlayer  continuity  conditions  between  layers  k  and  k-1  and  the  other  from  one 
interlayer  continuity  condition  between  k  and  k-1  and  a  layer  equilibrium  equation  involving  <|>k. 
Hence,  a  set  of  linear  simultaneous  equations  was  obtained  for  derivative  functions  (U^)' ,  (V^)' , 
(IJk)' ,  (Vk)',  etc.  in  terms  of  the  variables  IJk ,  Vk ,  Uk  ,  Vk  ,  etc.  These  functions  could  be 
easily  established  by  means  of  a  simple  simultaneous  equation  solver  on  the  computer.  The  IMSL 
library  subroutine  DVCPR ,  was  thus  utilized  for  obtaining  the  solution  for  Uk ,  Vk ,  etc. 

At  this  juncture,  it  is  important  to  note  how  the  boundary  conditions  at  the  free  edge  were 
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used  in  the  solution  procedure.  Again  the  weighted  displacemratt  parameter  <|>  was  eliminated  from 
the  appropriate  expressicms  ( the  expressions  for  Ny  and  N^y  in  this  case  ).  This  was  done  by 
obtaining  one  equation  (for  each  layer)  by  considering  simultaneously  Ny  and  N^y  at  y  =  +b 
and  again  at  y  =  -b.  The  final  boundary  equations  obtained  by  setting  y  =  ±  b  in  various 
expressions  such  as  Vy,  My  ,  Mjjy  ,  etc.  contained  the  weighted  displacement  parameters  and  the 
interhuninar  stress  components  specified  at  the  edges  y  =  ±  b. 

Finally,  there  were  (9N-2)  equations  obtained  by  the  layer  equilibrium  equations  and 
intersurface  continuity  and  top  and  bottom  surface  conditions.  To  these  were  added  6N  equations 
of  the  type  (2.1),  thus  adding  to  (15N-2)  equations.  The  unknown  parameters  were  given  as 
follows: 


tl  =  si  =  0  ;  P2N=  S2N  =  t2l^  =  0  Si  and  q.  . . .  2(N-1) 

U,  V,  W,  n ,  V ,  X  .  6N 

U,V,W,  ft,v,x  .  6N  (2.3) 

Pi  . (N-l)+  1 

Total  . (15N-2) 

In  using  the  library  subroutine,  (15N-2)  single  order  linear  differential  equations  were 
considered  for  solution.  The  number  of  boundary  conditions  are  listed  as  follows: 


Vyk  (±b)  =  0 

Myk  (±b)  =  0  . k=  1,2,3 . N 

M^^yk  (±b)  =  0 

S2k  (±b)  =0  ...  k=  1,2,3,..  .,N-1 

combination  of  equations  for  Nyk  and  N^yk  at  ±  b  _ k  =  1,2,3,. . .,  N 

equations  of  type  IJk'  -  Uk  =0  .  .  k=  1,2,3 . ,N 


Thus,  the  total  number  of  boundary  conditions  is  (15N-2). 


(2.4) 
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In  this  method  of  solution  however,  two  main  difficulties  were  encountered: 

(i)  The  reduction  of  the  numb^  of  equations  to  be  solved  by  elimination  made 
the  method  very  ctxi^licated  to  code, 

(ii)  The  Gauss- Jordan  computer  subroutine  that  reduced  the  Hrst  order  equations  to 
the  form  given  in  (2.2)  to  be  used  later  cm  in  the  differential  equation  solver 
DVCPR  encountered  problems  of  singularity. 

Thus  it  was  decided  not  to  eliminate  0 ,  but  instead  to  retain  the  original  set  of  equations 
because  the  effect  of  reducing  the  number  of  equations  to  be  solved  by  eliminating  0  was  not  found 
to  be  significant  enough  to  offset  the  additional  coding  involved.  The  cause  for  the  singularity 
arising  in  the  computer  subroutine  was  traced  to  be  due  to  the  coupling  of  the  displacement 
parameters  and  x*'  ht  the  layer  equilibrium  equation  6  of  reference  26.  This  necessitated  the 
introduction  of  another  variable 


Hk  =  wk  .  3^  =  Wk‘  -  x*"’  (2.5) 

This  substitution  had  an  additional  advantage  of  reducing  the  number  of  variables  to  be 
solved  for  in  each  layer  by  one.  After  carrying  out  the  above  substitution  and  judiciously 
manipulating  the  equations,  a  ncm-singular  system  of  equations  was  obtained  that  could  subsequently 
be  reduced  to  the  form  illustrated  by  equation  (2.2),  using  a  Gauss- Jordan  elimination  algorithm. 

If  the  interlaminar  equations 

S2*^  =  and  (2.6) 

tjk  =  tik+1 

were  incorporated  into  the  layer  equilibrium  equations  and  the  interlaminar  continuity  conditions, 
then  the  number  of  unknowns  in  each  layer  could  be  reduced  to  15,  viz.,  U,  V,  W,  G,  >)/,  x,  0, 
U,  V,  H,  12,  v ,  Pi ,  Si  and  ti  .  Additionally  there  were  three  extra  variables  P2^ ,  52^  and 
t2^  in  layer  N  (at  the  free  surface  of  the  laminate).  Since  the  laminate  was  unstressed  at  its  surface 
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pjN  «  sjN  =  tjN  =  0  (2.7) 

Furthermore,  in  the  first  layer,  due  to  the  symmetry  of  the  layer  about  the  midplane 

sil  =  tjl  =  0  (2.8) 

Hence  the  total  numbor  of  variables  to  be  solved  for  as  well  as  the  number  of  equations  became 
(15N-2). 


The  library  subroutine  DVCPR  required  that  there  be  the  same  number  of  boundary 
conditions  as  the  number  of  equations.  These  (lSN-2)  boundary  conditions  were  taken  as 

Vyk  (±  b)  =  Myk  (±  b)  =  (t  b )  =  0  . . .  k  =  1,2,. .  N 

Ny>^  (+b)  -  Nyk  (-b)  =  0 

N^yk  (+ b)  -  N^yk  (- b)  =  0  ..  k=  1,2,3,..  N-l 

Nyk(+b)  =*  N^yk(+b)  =  0 

S2k(±b)  =  0  ...  k  =  1,2,3,.  .N 

and  equations  of  type  Uk'  -  =0  . . .  k  =  1,2,3, . . .  N  (2.9) 

The  above  system  of  first  order  differential  equations  and  boundary  conditions  could  finally  be 
solved  by  the  previously  mentioned  library  subroutine. 

Even  in  this  method,  some  difficulties  were  encountered  due  to  the  fact  that  some  of  the 
layer  equilibrium  equations  and  the  interlaminar  continuity  conditions  had  to  be  differentiated  in  order 
to  achieve  a  form  similar  to  the  one  shown  by  equation  (2.1.2).  In  doing  so  the  original  problem  is 
weakened  since  the  order  of  some  of  the  equations  is  increased. 

Hence  in  order  to  preserve  the  nature  of  the  original  problem  and  to  avoid  singularity 
situation,  it  was  decided  to  construct  a  finite  difference  model  for  the  problem. 


14 


2.2  FINITE  DIFFERENCE  SOLUTION 


Since  the  free  edge  problem  had  been  reduced  to  essentially  a  one-dimensional  problem 
(in  the  width  dimension  y )  [  26],  the  finite  difference  model  was  based  on  formulas  for  the  first 
and  second  derivatives  derived  from  the  quadratic  Lagrange  polynomials  in  one  variable.  Since 
these  formulas  are  of  order  two,  the  enx)rs  vary  as  the  square  of  the  step  size. 

For  the  problem  on  hand,  variable  step  size  was  chosen  over  fixed  step  size.  This  was 
done  because  the  computations  could  then  be  made  more  efficient  by  making  the  step  size  at  regions 
of  high  stress  gradients  very  small  and  choosing  step  sizes  for  other  regions  acccxtiing  to  the  desired 
accuracy  of  the  results  desired.  Furthermore,  the  need  for  ficticious  nodes  adjacent  to  but  outside 
the  boundaries  are  eliminated.  This  was  accomplished  by  using  forward  and  backward  difference 
formulas  at  the  left  and  right  boundaries  respectively.  At  all  interior  nodes  central  difference 
formula  was  used. 

Finite  difference  solutions  to  similar  problems  have  been  given  while  studying  the  elastic 
response  of  involute  bodies  [27,28].  These  methods  are  direct  methods  used  to  solve  the  resulting 
set  of  finite  difference  equations.  Since  the  number  of  these  equations  was  quite  large,  the  limits  of 
available  computer  memory  was  attained  even  while  solving  relatively  small  size  problems.  In  order 
to  reduce  the  storage  burden  on  the  computers  and  to  make  the  calculations  more  efficient,  the 
iterative  Gauss-Siedel  method  was  used  in  the  current  solution  technique  to  solve  the  finite  difference 
equations.  Another  attractive  feature  of  using  the  iterative  method  was  that  by  reducing  the  memory 
requirement  of  this  procedure  the  maximum  problem  size  could  be  increased.  However,  this 
method  failed  to  yield  good  results. 

The  difficulty  seems  to  be  in  the  fact  that  the  way  the  original  problem  is  presented,  the 
number  of  boundary  conditions  is  less  than  the  number  of  field  equations  available  for  solution.  It 
has  been  found  [27,  28]  that  this  leads  to  inconsistencies  in  standard  numerical  methods.  Thus 
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some  of  the  field  equations  that  normally  are  satisfied  at  all  interior  nodes  may  be  enforced  at  the 
boundaries  so  that  they  supplement  the  boundary  conditions.  This  has  been  shown  in  the  above 
mentioned  references  to  yield  satisfactory  results. 
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3.  STIFFNESS  FINITE  ELEMENT  METHOD  FOR  GLOBAL-LOCAL 


MODEL 

The  global-local  variational  approach  has  already  been  numerically  solved  by  an 
exponential  series  solution  [25]  with  limited  success.  One  global-local  interface  had  been  considered 
in  that  procediffe.  In  the  present  section,  a  discretization  procedure  using  finite  element  method  of 
solving  the  laminate  problem  is  given  based  on  the  same  global-local  approach.  In  general,  instead 
of  considering  the  first  variation  of  the  energy  functional  representing  the  laminate,  the  stresses  and 
displacements  are  written  in  terms  of  the  nodal  degrees  of  freedom  of  an  element  of  chosen  shape. 
The  solution  to  the  problem  is  then  sought  by  making  the  functional  stationary  with  respect  to  each  of 
the  nodal  degrees  of  freedom.  At  first,  a  current  literature  survey  was  conducted  pertaining  to  the 
finite  element  models  based  on  mixed  variational  principles  (Reissner's  principle,  in  this  case).  The 
purpose  of  this  literature  review  was  to  establish  the  guidelines  for  the  application  of  the  present 
finite  element  technique  to  the  global-local  model  of  the  laminate.  The  variational  functional  for  the 
entire  laminate  consists  of  two  distinct  energy  functionals  -  one,  a  total  potential  energy  functional 
characterizing  the  global  region  and  the  other,  a  Reissner  variational  functional  for  a  predetermined 
local  region  of  interest  (figure  1).  The  finite  element  discretization  based  on  these  variational 
principles  is  individually  well  documented  [38,  39, 40, 41].  The  potential  energy  principle  is  very 
commonly  used  as  the  basis  for  a  displacement  -based  compatible  finite  element  model.  The 
Reissner  variational  principle  is  intrinsically  multifield  in  which  the  primary  field  variables 
characterizing  the  element  behavior  are  both  stresses  and  displacements. 

Following  the  theoretical  development  for  the  global-local  model,  the  variational 
functional  for  the  entire  laminate  is  written  to  represent  both  i^ie  global  region  and  N  layers  of  the 
local  region  as  follows: 

N 

^  =  /  Wp  (Uj ,  ejj )  dv  +  S  /  [  1/2  ajj  (uy  +  uj^j)  -  WcCOy  ,eij)]lf  dv^ 

Vg  k=l  Vk 

-  JtiUids  (3  1) 
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In  the  above  expression,  Wp  (u{ ,  ejj)  is  the  potential  energy  of  the  global  region  and 
Wj.  (ojj  ,  ejj )  is  the  canplemcntary  energy  corresponding  to  the  local  region.  The  terms  Vg  and  Vj^ 
refer  to  the  volume  domains  of  the  global  and  local  regions  respectively.  It  should  be  noted  that  Vj^ 
in  particular  refers  to  a  local  layer  k  in  the  local  domain  which  is  assumed  to  have  N  local  layers. 

In  the  abbreviated  form,  the  above  expression  can  be  written  as 

7C=  7Cg+7C/  +  JCs  (3.2) 

where  K  g ,  referring  to  the  first  integral  of  (3.1),  is  the  potential  energy  for  the  global  region,  and 
7C/  is  the  Reissner  variational  functional  for  N  layers  in  the  local  region.  7C  /  is  the  summation 
term  (k  =  1  to  N)  in  the  equation  (3.1).  Finally,  Kg,  referring  to  the  last  term  of  (3.1)  corresponds 
to  the  potential  energy  of  the  prescribed  surface  tractions. 

3. 1  ELEMENT  MODELING  FOR  THE  GLOBAL  REGION 

By  defining  the  constitutive  laws  as 

Oj  =  Cij  (  ejo  +  z  KjO  -  Cj )  (  i,j  =  1,2.3,6  ) 

Oj  =  Cij  (  ejo  +  z  Kj  +  z2/2  PjO  )  (  ij  =  4,5  )  (3.3) 

the  potential  energy  for  the  global  region  is  written  as 

7Cg  =  ^gE  (3.4) 

In  the  above  expressions,  e®  and  e  are  midsurface  and  expansional  strains  respectively, 
and  K  and  P  are  curvature  parameters  characterizing  strains  at  any  distance  z  from  the  surface. 
These  through-the-thickness  variations  of  strains  follow  those  for  displacements  on  the  global  region 
(i.e.,  a  linear  variation  for  inplane  displacements  and  a  quadratic  variation  for  the  transverse 


displacement).  TCgg  and  Kg^  in  (3.4)  refer  to  the  volume  integrals  over  volume  domain  Vg  of 
th  global  region  in  terms  of  the  mechanical  and  expansional  strains  respectively.  Because  of  the 
polynomial  variation  with  respect  to  z ,  of  the  displacements  (and  hence  of  strains),  it  is  possible  to 
write  each  of  the  integrals  TCg^  and  TCg^  as  a  sum  of  several  integrals  with  various  powers  of  z 
appearing  as  coefficients.  In  other  words,  it  is  possible  to  carry  out  the  integration,  with  respect  to 
z,  of  the  stiffness  terms  Qj  in  conjunction  with  various  powers  of  z  to  obtain  the  moduli 
corresponding  to  the  entire  global  domain.  These  moduli  are  defined  as 

Wl 

(  Ajj  ,  Bjj,  Djj  ,  Fjj ,  Hjj  )  =  J  (  1  ,  z  ,  z2  ,  z3  ,  z4  )  Qj  dz  (3. .‘5) 

-H/2 

where  H,  appearing  in  the  integration  limits,  is  the  total  thickness  of  global  domain  of  layers.  This 
process  reduces  the  volume  integral  to  an  area  integral  in  terms  of  the  x  and  y  coordinates.  For 
this  purpose,  the  various  components  of  strains  and  curvatures  are  written  in  vector  form  as  follows; 

{e9}T  =  {  d®  62°  630  egOjT 

{  eo}T  =  {  64°  ejOlT 
{  K  }T  =  {  Ki  K2  K3  Kg  }T 

{  K  }T  =  {  K4  K5  )T 

{  P  I'r  =  {  P4  Ps 

{  e  }T  =  {  ei  0203  eg  )T  (3-6) 


can  be  written  as 


Then  the  integral  corresponding  to  the  mechanical  strains  in  the  global  region 


7Cge=  1/2  J(  TC^gg  +z  7C2g£  +  z‘2n^g^+  T>  K  ^g^  )  dv 


(3.7) 


where  for  example 
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JCV  =  {eoF[Cii](eo}  +  }T[C22]{eo}  +  2{eo)T[C2i]{e»}  (3.8a) 


7c2ge  =  2{  K}T[C„]{eP}  +  2{  K}T[C22]{eO}  +  2{  K}T[C2,]{eO)  + 

2{  eo)T[C2i](ic}  (3.8b) 


and  so  on. 


The  material  stiffness  matrix  in  the  above  expressions  are  defined  as 


Cii 

Ci2 

Cl3 

C16 

[c,l]  = 

Cl2 

C22 

C23 

C26 

Ci3 

C23 

C33 

C36 

_Ci6 

^26 

^36 

^66 

[^]  = 

— 

Ci4 

Qs 

Qs 

C55_ 

[^2,]  = 

Ci4 

C24 

C34 

C46 

Cl5 

C25 

C35 

C56 

In  a  similar  way,  the  expression  for  TCgg  can  also  be  written  as 


TCge 

where  for  example, 

It'gj  =  (e<>  )T[C,i](e|  +  (e»  )T(C2,]{e) 

icV  =  {K  )T(Ciil{e)  +  (  K)T[C2,l{e) 


ge 


+  z^n 


ge 


)  dv 


(3.9) 


(3.9a) 

(3.9b) 
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For  the  finite  elemoit  description,  the  vector  of  variables  for  the  global  region  is 


{  u  =  {  u ,  V ,  u* ,  v*.  w  ,  w*  ,  w  (3.10) 

where  u  ,  v  ,  etc.  are  the  weighted  displacement  parameters  that  are  related  to  the  midsurface 
displacements  and  curvature  terms  [24]  by  the  following  relations: 

uo  =  1/2  u 

vO  =  1/2  V  etc.  (3.11) 

and  by  definititMi, 

H/2 

( w ,  w*,  w )  =  I  w  (  1 ,  2z/H2  ,  4z/h2  )  2/H  dz 
-H/2 

It  is  now  possible  to  write  each  of  the  vector  quantities  in  (3.6)  in  terms  of  the  vector  of 
variables  {  u  }.  For  example. 


{e«}  =  [Li1[L2]{u} 

(3.11a) 

(e»)  =  [Li][L2l{u) 

(3.11b) 

(K)  =  [L3][L4]{u} 

(3.11c) 

{  K  }  =  {  L3)  w* 

(3.1  Id) 

(P)  =  {L5}[L6]{u) 

(3.11c) 

In  the  above  expressions,  there  are  two  other  vectors  {  u  )  and  {  u  }  on  the  right  hand 
side  ( in  addition  to  a  single  weighted  displacement  parameter  w* )  which  arc  but  parts  of  the  vector 
of  complete  variables.  That  is. 


{  u  I’l'  =  {  u*  V*  w  w  }T  (3.12a) 

(  u  =  {  w  w  (3.12b) 

All  the  other  matrices  spearing  in  the  set  of  equations  (3.1 1)  are  given  below. 
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[L,] 


5/dx  0  0 

0  a/dy  0 

0  0  1 

a/dy  a/dx  0 


(3.13a) 


[L2] 


[L2] 


[LJ 


0 

0 

0 


1/2  0  0  0  0 

0  1/2  0  0  0 

0  0  0  0  0 

=  a/dy  0  1 

d/dx  1  0 

0  0  0  9/8 

0  3/H  0  0 

0  0  3/H  0 

[L3]  =  [L,] 


0  0 
0  0 
3/H  0 


0  -15/8 
0  0 

0  0 


(3.13b) 


(3.13c) 


(3.13d) 


(3.l3e) 


[L4] 


3/H  0  0  0 

=  0  3/H  0  0 

0  0  -15/H2  45/H^ 

{  L3  )  =  3/H  {  d/dy  d/dx 

{  L5  }  =  { d/dy  d/dx  jT' 


(3.13f) 

(3.l3g) 

(3.13h) 
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-15/H2 

0 


0 

45/H2 


(3.13i) 


[m  = 


3.2  DISCRETIZATION  PROCEDURE  FOR  THE  GLOBAL  REGION 


The  nature  of  the  interpolation  functions,  which  describe  the  mapping  of  the  variables 
chosen  for  the  element,  is  derived  from  the  governing  equations  for  the  global  domain.  The 
variables  are  given  in  (3.10).  From  the  stress  -  strain  relations  (3.3),  it  is  possible  to  obtain  the 
stress  -  displacement  relations  by  using  the  kinematic  relatioships  [  25  ]  given  below: 


EjO  =  du®  /  3x 
^  =  dvo  /  3y 


£3°  =  Vz 


(3.14) 


=  ¥z  3wO/dy  etc. 


Then,  it  can  be  seen  that  for  the  stresses  to  be  linear  in  the  spatial  coordinates  x  and  y, 
the  displacements  and  the  curvature  terms  vP,  vO,  w®,  \|/;j ,  Vy  ,  and  <J>  have  to  be  quadratic 
in  X  and  y.  Consequently,  from  the  relationships  (3.1 1) ,  the  weighted  displacement  parameters 
(3.10)  should  vary  quadratically  with  respect  to  x  and  y. 


In  the  present  formulation,  triangular  elements  are  chosen  for  discretization  purposes.  It 
has  been  established  [38,39]  that  it  is  convenient  to  work  with  area  coordinates  or  natural  coordinates 
when  dealing  with  triangular  elements.  If  ,  ^2  and  ^3  are  the  area  coordinates,  then  the 
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relationships  between  these  and  the  cartesian  coordinates  are  given  by  ( figure  2.)  the  following  set 
of  equations.  The  last  of  the  following  expr^sions  is  an  identity  associated  with  the  area 
coordinates.  This  identity  is  used  to  obtain  a  unique  relationship  between  the  two  sets  of 
coordinates: 


X  =  xi  Cl  +  *2  C2  +  X3  C3 
y  =  YiCi  +  y2C2  +  ysCs  (3-i5a) 

1  =  Cl  +  C2  +  C3 


Following  the  above  equations,  it  should  be  noted  that  the  edge  1-3  of  the  triangle  is 
deHned  by  only  two  area  coordinates  Ci  C2  <  while  the  third  area  coordinate  C3  =  0. 
Similarly,  the  boundary  3-5  has  Cl  =  0  ®nd  on  boundary  5-1 ,  the  coordinate  C2  =  0*  ^so, 
the  numerical  integration  of  functions  on  the  rectangular  boundary  1-3  is  carried  with  respect  to  Ci 
and  C2  oi^y*  Similar  reasoning  follows  for  the  otiier  two  rectangular  boundaries  and  the  integraticm 
is  canied  out  with  respect  to  aU  the  three  coordinates  for  the  upper  and  lower  triangular  boundaries 
of  the  element 


By  inversion  of  equations  (3.15a),  one  can  get 


'Ci' 

y23 

X32 

*2  y3 

-  ^^3  y2 

C2 

►  =  1/2A 

y3i 

*13 

X3  yi 

*  xi  ys 

1^^ 

yi2 

*21 

*1  y2 

-  *2  yi 

(3.15b) 


and 


In  the  above  equations,  x^,  yj  (i  =  1,2,3)  are  the  cartesian  coordinates  of  the  node  i 


y23  =  y2  -  y3 
X32  = 

^^21  = 


X3  -  X2 

X2  -  Xj  etc.  and 
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2A  =  (xi-  X3)(y2  -  ys)  -  (X2  -  X3)(yi  -  y3) 


For  the  puipose  of  establishing  rules  for  differentiation  and  integration,  only  and  ^2 
are  considered  as  independent  variables.  Thus  d /dCi  implies  that  ^2  is  held  constant  and  d/d^2 
implies  that  is  held  constant  The  differentiation  rule  is  given  by 


(3.16a) 


where  [  J  ]  is  the  Jacobian  matrix.  Its  determinant  is  equal  to  2A.  In  other  words 


[J]  “ 


*13 

yi3 

*23 

y23 

Thus  the  differential  area  dxdy  is  replaced  in  all  integrations  by 


dxdy  =  IJI  d^i  dl^2 


(3.16b) 


(3.16c) 


where  IJI  symbolizes  the  determinant  of  the  Jacobian  Matrix  [J]. 


In  the  final  analysis,  it  is  necessary  to  maintain  continuity  of  displacements  u  ,  v  and  w 
between  the  layers  of  the  laminate.  For  this  puipose,  we  choose  as  the  nodal  degrees  of  freedom  the 
displacements  at  the  top  and  the  bottom  of  the  layer.  (Denoting  by  the  subscripts  t  and  b  the  top 
and  the  bottom  displacements  respectively  (such  as  u^  and  U|,  etc.),  it  is  possible  to  use  the 
relationships  such  as,  for  example. 


u  =  1/2  u  +  3z/H  u* 
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and  write 


Ut  =  1/2  u  +3/2u* 

Ub=  1/2 Q  -3/2u* 

vt  =  l/2v  -3/2v*  (3.17a) 

Vb=  1/2  V  -3/2v* 

Wj  =  -3/4  w  +3/2  w*  +  15/4  w 
Wj,  =  -3/4  w  -  3/2  w*  +15/4  w 

Another  relationship  fw  w  is  added  so  that  the  seven  unknown  displacements  are  related 
to  the  seven  weighted  displacement  parameters.  This  last  relationship  refers  to  the  middle  surface 
displacement  Wq  .  That  is, 


Wq  =  9/8  w  - 15/8  w  (3.i7b) 

The  displacements  u^,  uj,,  Vj,  vj,,  Wj,  wi,,  Wq  are  graphically  shown  for  different 
nodes  of  a  layer  in  Figure  2. 

In  the  matrix  form,  the  equations  (3.17)  can  be  written  as 

(utb)  =  [Luul  {«)  (3.18) 

where  {  Utb)  =  {  Uj,  Ub,  Vj,  Vb,  Wj,  Wb,  Wg  )T  and  the  elements  of  the  matrix  [  Luq]  are 
given  by 
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1/2  0  3/2  0 

0  1/2  0  3/2 

1/2  0  -3/2  0 


0 


0 


0 


0  0  0 

0  0  0 

[  Luq  ]  =  0  1/2  0  -3/2  0  0  0 

0  0  0  0  -3/4  3/2  15/4 

0  0  0  0  -3/4  -3/2  15/4 

0  0  0  0  9/8  0  -15/8 

The  inverse  relationships  of  (3.18)  can  be  writtoi  as 

{u}  =  [Luq]-Mu*}  (3.19) 


As  per  earlier  discussion,  the  inteipolation  representation  for  the  displacements  {  Ufi,}  is 
written  as,  fw  example, 

Uj  *  Ujl  Nj  +  N2  +  0(3  N3  +  Ui**  N4  +  UjS  N5  +  Uj®  Ng  (3.20) 

Similar  quadratic  representations  are  written  for  ug ,  Vj ,  vg  ,  wg ,  Wj,  also.  In  the 
above  expression  (3.20),  Nj  ( i=l,6 )  represent  the  shape  functions  and  Uj*  and  ug*  etc.  refer  to 
the  corresponding  quantities  at  the  node  i.  The  same  shape  functions  are  used  to  charact^e  all  the 
displacements.  These  shape  functions  are  given  by  the  following  expressions  [  38  ]: 

Ni  =  Ci  (2Ci  -  1)  (i=  1.2.3) 

N4  =  4  C1C2 

N5  =  4  C2C3 

N6=4CiC3  (3-21) 

It  is  now  possible  to  use  equations  (3.18)  through  (3.21)  in  conjunction  with  equations 
(3.11)  to  obtain  the  various  stiffness  matrix  terms  for  the  global  region.  As  an  illustration,  the 
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following  development  is  given. 


FromO.lla),  {  e® }  =  [Lj]  [LjKu) 

With (3.19),  {  e®}  =  [Lil  [L2I  [L^rMu*} 

With  the  set  of  interpolation  functions  (following  the  type  given  in  (3.20))  given  as 

(u*)  =  [NsllUeN)  (3.22) 

in  which  the  vector  {  }  includes  all  the  nodal  degrees  of  freedom  for  the  element.  Then, 

{  e®}  =  [Li][L2][Lj,fi]-MNs]{UeN} 

=  [LiiNKUeN)  (3.23) 

Hence  the  first  expression  on  the  right  hand  side  of  equation  (3.8a)  is  written  as 

1  e“lTtCiil(  E”)  -  (ll«N)T(L,jN)TtC„][L,2N](UeN)  (3.24) 

The  corresponding  stiffness  matrix  is 

J  [LijN]  [Cu]  [LijN]  dv 

Vg 

In  the  above  expression,  the  integration  with  respect  to  z  can  be  done  independently. 
The  remaining  terms  in  the  integrand  are  now  functions  of  x  and  y  only.  These  area  integrals  are 
numerically  evaluated  by  using  the  Gauss  quadrature  formulas.  Similar  expressions  as  (3.24)  are 
written  for  various  quantities  given  in  equations  (3.7)  and  (3.8)  and  element  stiffness  matrices  are 
individually  written  for  each  expression  of  these  equations.  In  each  case,  the  integration  with 
respect  to  z  is  appropriately  done  by  considering  the  various  powers  of  z  appearing  in  (3.7). 
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3.3  ELEMENT  MODELING  FOR  THE  LOCAL  REGION 


The  functional  characterizing  the  behavior  of  the  entire  laminate  also  consists  of  a 
Reissner's  functional  that  represents  the  local  region  of  the  laminate.  For  a  local  region 
consisting  of  N  layers,  this  functional  is  given  as 

N 

JCr  =  S  JC/*'  -  J  X  U  ds  (3.25) 

k=l  S<y 

where  7C/  ^  is  the  volume  integral  expression  given  by 

JC/k  =  /  [1/2  ay(uy  +  uj4)  -  U*(a^,eij)]k  dvk  (3.26) 

Here,  U*  ( Ojj ,  ejj )  is  the  complementary  energy  density  and  V/  ^  refers  to  the  volume  enclosed 
by  the  layer  in  the  local  region.  The  volume  integral  expression  is  written  as  (without  the 
superscript  k  for  convenience), 

TC/  =  J  [  a  (  e  -  e )  -  U*  ]  dv  (3.27) 

V/ 

Writing  e  =  ( S  o  +  e )  and  the  appropriate  expression  for  U*,  the  equation  (3.27)  can 
be  rewritten  as 


Ki  =  1/2  J  aT  S  0  dv  +  J  0T  e  dv  (3.28) 

V/  V, 

where  S  represents  the  compliance  matrix.  It  should  be  noted  that  o,  e  and  e  are  vector  quantities. 
For  the  local  region,  the  in-plane  stress  components  are  assumed  to  be  linear  functions  of  z  and  the 
remaining  stresses  are  either  quadratic  or  cubic  in  z  as  a  consequence  of  substitution  in  the  layer 
equilibrium  equations.  These  stresses  are  functions  of  the  interlaminar  stress  components  pj ,  P2 , 
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t| ,  t2 ,  sj  and  S2  [24].  These  latter  terms,  therefore  appear  as  the  nodal  degrees  of  freedom  for 
the  layer  in  the  local  region  in  addition  to  the  displacement  degrees  of  freedom  as  in  the  global 
region. 


As  in  the  case  of  strain  components  for  the  global  region,  the  stress  components  here  are 
split  up  into  two  vectors  for  mathematical  convenience.  These  are 

{  <yi  )T  -  {  Oj  (J2  03  06  and  {  o2  }T  =  {  04  O5  }f  (3.29) 


By  using  relations  (3.25)  of  reference  24,  these  vectors  can  be  written  in  terms  of  the 
force  and  moment  components  and  the  interlaminar  stress  components.  These  are 


[oM=[Mi*1{Nr}  and  (o2}  =  [M2*l{V]  (3.30) 

where 


{  Nr  =  { Njj/h^  Ny/h2  Nj/h^  Mjj/h^  My/h^  M^jy/h^  pj  P2 )  ^  » 

{V}T  =  [V^/h  Vy/h  Si  S2  ti  t2}  and 


[Ml*] 


10  0  0 

0  10  0 

0  0  Fi  0 

0  0  0  1 


12z/h  0  0 

0  llzJYi  0 
0  0  F2 

0  0  0 


0  0  0 

0  0  0 

0  F3  F4 

12z/h  0  0 


1  [  M2*  ] 


0 

Fl 


Fi  F5  F6  0  0 

0  0  0  F5  Fg 


In  the  above  expressions,  quantities  Vx ,  Vy ,  Nj, ,  M^  are  deHned  in  reference  24,  h  is 
the  thickness  of  the  layer  under  consideration  and  Fi  through  Fg  are  defined  as  follows: 
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Fi  = 

3/2  - 

6z2/h2 

II 

30z/h 

-  120z2/h2 

II 

-1/4  + 

3z/2h 

+  3z2/h2 

-  10z3^3 

F4  = 

-1/4  - 

3z/2h 

+  3z2/h2 

+  10  z3/h3 

F5  = 

-1/4  - 

z/h  + 

3z2/h2 

and 

F6  = 

-1/4  + 

z/h  + 

3z2/h2 

Corresponding  to  the  stress  components  (3.21),  the  compliance  matrix  for  the  layer  is 
split  up  as  follows: 


[s.l]  = 

[§22] 

[S21]  = 

Again  as  before, 

7C/  =  ^  -f 


Sii 

S12 

Sl3 

S16 

S12 

S22 

S23 

S26 

S13 

S23 

S26 

S36 

S16 

S26 

S36 

^66 

S44 

S45 

S45 

S55 

Si4 

S24 

S34 

S46 

Sl5 

S25 

S35 

S56 

2  +  7C/3  +  7C/4  + 


(3.31) 


where 

7C/1  =  1/2  j  {  )T[Sn]  {  )  dv  and 
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7C,2  =  1/2  J  {a2}T[S2i]{ol}  dv  and 


so  on. 


Upon  substitution  of  expressions  (3.30)  for  {  }  and  {  )  into  equations  of  type 


(3.32),  one  gets,  for  example. 


1  =  1/2  J  {NR)T[Mi*F[Sn]  [Ml*]  {Nr}  dv 


It  should  be  noted  that  in  the  above  expression,  the  z-terms  appear  only  in  the  matrix 
[M 1  *] .  As  such,  the  above  expression  can  be  rewritten  as 


.1  =  1/2  J  {Nr}T[S,i*]  (Nr)  dA 


(3.33) 


where,  now. 


[Sii*]  =  J  [Mi*]T[S,i]  [Ml*]  dz 


The  other  expressions  in  (3.32)  can  be  similarly  written  as,  for  example , 


=  1/2  I  {  V}T  [S21*]  [Nr]  dA  etc. 


(3.34) 


Thus  it  is  possible  to  write  each  integral  expression  in  equation  (3.31)  in  terms  of  the 
force  and  moment  components  as  well  as  the  interlaminar  stress  components. 


3.4  DISCRETIZATION  PROCEDURE  FOR  LOCAL  REGION 


It  is  now  possible  to  invoke  the  layer  constitutive  relations  to  write  (3.34)  in  terms  of 
the  displacement  functions  and  the  interlaminar  stress  components.  These  can  be  symbolically 


wntten  as 
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{Nr)  =  [RiKui)  +  [RsKel 


(3.35a) 

(3.35b) 


(V)  =  [RjKuj) 

where 

{ Uj)  =  { u  V  u*  V*  w  w*  w  P]  P2  } 

{ U2}  =  { u*  V*  w  w  Si  S2  ti  t2 } 

The  matrices  [  Ri) ,  [  R2]  and  [  R3]  express  the  relations  of  the  force  and  moment 
components  with  the  above  mentioned  vectors  of  the  weighted  displacement  parameters  and  the 
interlaminar  stress  components.  Substitution  of  expressions  (3.35)  into  expressions  of  type  (3.33) 
and  (3.34)  results  in  the  variational  functional  1C/  of  (3.31)  to  be  written  as 

It,  =  1/2  J[(tt,)TtR„)(u,  )+  {u,  )T  [R,3l  le)  + 

{«)T[R„JT{„,)  +  {e)T(R33]  (e)  + 

(U2FtR22)  (U2)  +  (U|)T(R,2](U2)  + 
(U2)TtRl2l’'(Ui)+  (Uj)’’ [R23I  (e)  + 

(e)T  [R23lT{„2)  ]  dA  (3.36) 


The  various  matrices  [Rnl,  [RiJ  etc.  are  defined  as  follows: 

[Rill  =  [Ril  [Sii*]  [Ri] 

[R12]  =  [Ri]TtS2i*]T[R2] 

[R23]  =  [R2]’^[S2i*]  [R3]  and  so  on. 


Again,  as  in  the  case  of  global  region,  it  is  necessary  to  maintain  displacement  and  stress 
continuity  between  any  two  layers  of  the  local  region.  Thus,  we  revert  to  actual  set  of  degrees  of 
freedom  which  are  now  the  actual  displacements  at  the  top  and  bottom  of  each  layer  ( u^ ,  uj,  etc.) 
along  with  the  interlaminar  stress  components. 
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The  relationships  between  the  actual  displacements  and  the  weighted  displacement 
parameters  are  similar  to  equations  (3.17)  [24]  and  the  steps  fm*  further  development  of  element 
stiffness  matrices  for  the  local  region  follow  those  for  the  global  region.  The  major  differences 
between  these  two  developments  are  that: 

(i)  in  die  matrix  expression  in  equation  (3.18),  there  are  additional  terms  refaring  to  the 
interlaminar  stress  components, 

(ii)  while  the  interpolation  functions  for  the  displacements  follow  the  same  patt^  of 
quadratic  variation,  the  new  degrees  of  freedom  for  the  local  region  referring  to  the  interlaminar 
stress  components  are  linearly  interpolated.  That  is. 


Pi  = 

Pl‘N, 

+  Pi2N2  +  Pi3N3 

P2  = 

P2'Ni 

+  +  P2^N3 

tl  = 

t,>N, 

+  ti2N2  +  ti3N3 

*2  = 

IjlNi 

+  t22N2  +  t22N3 

Si  = 

S.IN, 

+  Si2N2  +  S,3N3 

S2  = 

SjlN, 

+  S22N2  +  S23N3 

(3.37) 


These  expressions  are  incorporated  into  the  equation  (3.36)  and  the  variation  of  7C{ 
with  respect  to  the  chosen  degrees  of  freedom  generates  the  stiffness  matrix  for  the  element  of  the 
local  layer  under  consideration. 


3.5  surfac:e  integrals  involving  local  and  global  regions 


The  surface  integrals  to  be  considered  in  the  present  variational  formulation  for  the 
global-local  model  are 

N 

JfudS+JxudS  +  01^+* )  dSij  +  |  ( X/  +  Xg )  u  dS 

So/I  S^gi  kal  Sfjii  SQg2 

(3.38) 
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In  the  above  expression,  ^o/l  and  ^agl  represent  in  general  the  external  boundaries  of 
the  local  and  global  regions,  req>ectively  where  the  tractions  are  prescribed.  These  include  the  top 
and  bottom  surfaces  of  the  laminate  (where  an  external  loading  may  or  may  not  be  prescribed).  In 
(3.38),  the  additional  subscript  2  represents  the  interlaminar  surfaces  between  the  local  layers  and 
the  commcM)  surface  between  the  global  and  local  regions.  Since  there  are  no  prescribed  tractions  on 
^0/2  and  Sag2  (the  continuity  of  stresses  and  displacements  being  taken  care  of  by  common  nodal 
parameters),  the  last  two  integral  terms  of  (3.38)  are  zero.  Any  of  the  two  remaining  integrals  can 
be  generally  written  as 


I  t  u  dS  (3.39) 

So 

where  {t  }  is  the  traction  vector  of  boundary  tractions  on  the  boundary  S(f  and  {  u  }  is  the  vector 
of  corresponding  displacement  components.  These  displacement  compcments  can  be  written  in  terms 
of  the  appropriate  nodal  degrees  of  freedom  pertaining  to  the  boundary  under  consideration.  The 
equation  (3.39)  can  include  a  general  type  of  loading  that  can  be  represented  by  the  vector  {f  } .  The 
traction  components  in  the  surface  integral  expression  need  to  be  specified  on  every  boundary  of  the 
element.  By  a  knowledge  of  the  state  of  stress  at  any  location,  the  traction  components  on 
boundary  surfaces  are  readily  obtained  by  the  following  general  relatioships. 


''xy 

“^XZ 

\ 

''xy 

Oy 

'^yz 

)  V 

(3.40) 

( "v 

In  the  above  equations  n^  ,  ny ,  n^  are  the  direction  cosines  of  the  normal  to  the  boundary 
surface  under  consideration. 

In  the  present  mathematical  formulation,  a  uniformly  distributed  load  is  considered  for 
checking  the  validity  of  the  present  approach.  In  such  a  case  of  transverse  loading,  only  the  z- 
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component  of  the  traction  vector  is  non-zero  and  is  equal  to  the  magnitude  of  the  applied  load  (n2  =  1 
and  n^  =  Uy  =  0).  The  surface  integral  (3.39)  then  beccmies 

q  J(N}i,T{u}b  dS  (3.41) 

Sa 

where  (  N  }b  refers  to  the  vector  of  shape  functions  associated  with  the  dispacement  w  and  {  u  )b 
is  the  vector  of  the  degrees  of  freedom  pertuning  to  the  boundary.  This  expression  now  involves 
integration  with  respect  to  x  and  y  of  the  boundary.  As  in  the  previous  cases  this  integration  is 
written  in  terms  of  the  area  coordinates  and  the  numerical  integration  is  done  by  Gaussian  quadrature 
constants  as  described  in  the  next  section.  This  results  in  a  consistent  load  vector  as  the  right  hand 
side  of  the  final  equations  to  be  solved. 

3.6  NUMERICAL  SOLUTION  PROCEDURES 

By  considering  the  mathematical  developments  of  discretization  as  in  the  previous 
sections,  a  numerical  solution  scheme  was  evolved  and  the  corresponding  computer  algorithm 
written.  This  corresponds  to  the  various  stages  such  as  evaluation  of  volume  integrals  for  global 
and  local  regions,  assembly  procedures  at  the  element  level  and  later  at  the  structure  level  for  the 
triangular  mesh  pattern  chosen  and  the  appropriate  inclusion  of  boundary  conditions  and  finally 
solution  of  algebraic  equations.  As  mentioned  before,  the  initial  cases  considered  for  the  numerical 
work  pertain  to  the  checking  of  the  validity  of  the  present  model.  As  such,  numerical  problems  for 
which  solutions  are  established  are  chosen  for  this  purpose. 

The  final  from  of  variational  functional  governing  the  behavior  of  the  laminate  can  be 

written  as 


TT  =  1/2  K  f  -  fTQ 


(3.42) 


In  this  expression,  the  stiffness  matrix  K  is  obtained  by  assembling  the  corresponding 
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stifChess  matrices  for  all  the  elements.  Thus  K  can  be  referred  to  as  the  structure  stiffness  matrix. 
Similarly,  the  consistent  load  vector  Q  also  corresponds  to  the  structure,  f  in  the  above 
expression  is  the  vector  of  die  nodal  degrees  of  fteetkim  for  the  structure.  The  variation  of  the  above 
functional  with  respect  to  each  of  the  nodal  degrees  of  freedom  gives  rise  to  a  set  of  linear, 
simultaneous  equations.  It  is  then  necessary  to  apply  the  boundary  conditions  applicable  to  the 
problem  under  consideration.  This  process  reduces  the  number  of  simultaneous  equations  to  be 
solved  finally. 

The  various  analytical  developments  given  in  the  previous  sections  have  been 
transformed  into  computer  codes  in  order  to  solve  numerical  problems  pertaining  to  laminate  stress 
analysis.  The  computer  program  written  for  this  purpose  covers  the  following  topics.  The 
cmresponding  subroutines  perfonning  these  qierations  are  also  given. 

(1)  Input  element  description,  nodal  coordinates  and  properties  of  individual  layers  ( 
thickness,  fiber  orientation,  material  constants  etc.).  The  nodal  number  generation  and  the 
corresponding  coordinate  evaluation  are  done  by  the  subroutine  ELGEN  after  the  dimensions  of  the 
plate  and  the  mesh/grid  patterns  are  read. 

(2)  Calculation  of  global  modulii  for  the  layers  comprising  of  the  global  region.  This  is 
achieved  by  the  subroutine  GLOBQ,  the  output  of  which  are  the  quantities  of  the  layer  moduli 
integrated  appropriately  through  the  thicknesses  of  the  layo^. 

(3)  Calculation  of  shape  function  parameters  and  Jacobian  of  transformation  which  are 
applicable  to  both  global  and  local  regions. 

(4)  Evaluation  of  the  various  stiffness  matrix  elements  for  global  region  using  the  global 
moduli.  This  is  done  in  a  subroutine  VOGLOB  which  essentially  converts  volume  integral  terms  to 
area  integral  terms.  The  numerical  integration  is  effected  by  the  5-point  Gauss-Radau  formula  for 
which  the  constants  are  read  as  input  in  the  beginning. 
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(5)  Assembly  of  the  matrix  elements  at  the  element  level  for  the  global  region  using 
subroutine  ELSTG. 

(6)  Evaluation  of  compliance  elements  for  a  local  layer  based  on  the  individual  layer 
properties  for  that  layer.  This  evaluation  is  done  in  the  main  program.  Calculation  of  stiffness 
matrix  elements  for  a  local  layer  following  the  analytical  descriptions  is  also  done  here. 

(7)  Assembly  of  stiffness  matrix  elements  at  the  element  level  for  the  local  layer  under 
consideration  (subroutine  ELSTL). 

(8)  Repetition  of  the  calculations  for  all  the  local  layers  so  that  the  stiffness  matrix 
for  the  element  now  includes  all  the  local  layers  as  well  as  the  global  region.  Repetition  of 
calculations  of  the  above  steps  to  cover  all  the  elemoits  of  the  structure. 

(9)  Calculation  of  the  consistent  load  vector  corresponding  to  the  element  under 
consideration  (subroutine  LOAD). 

(10)  Assembly  of  element  stiffness  matrix  and  load  vector  elements  to  obtain  the 
corresponding  structure  stiffness  matrix  and  load  vectm*  using  the  subrouitne  ASSMB. 

(11)  The  boundary  conditions  read  in  step  (1)  are  now  invoked  (subrouime  BCOND) 
and  the  necessary  rows  and  columns  are  set  to  zero  corresponding  to  zero  boundary  conditions. 

(12)  Solution  of  the  resulting  set  of  simultaneous  equations  to  obtain  the  required  nodal 
degrees  of  freedom.  This  is  done  by  means  of  library  subrouitne  LEQTIP.  The  results  -  both 
stresses  and  displacements  -  are  listed. 

These  steps  are  shown  in  the  form  of  a  flow  chart  in  Figure  4. 


38 


As  described  in  the  previous  sections,  area  integrals  involving  cartesian  coordinates  are 
encountered  during  the  present  discretization  process.  Since  the  area  coordinates  ( or  the  natural 
coordinates )  are  chosen  for  die  numerical  work,  the  computations  are  also  dcme  with  the  latter  as  the 
variables.  Specifically,  a  typical  integral  in  cartesian  coordinates  is  transformed  to  one  in  area 
coordinates  as  follows. 

I  =  JJ  F  dxdy  =  JJ  F  detlJ]  d^i  dC2 


The  numerical  evaluation  of  the  above  integral  is  carried  out  by  means  of  the 
Gauss-Radau  integration  constants  involving  5  points  [38].  In  other  words,  the  integral  is 
numerically  written  as 


where 


5  5 

I  =  EX  W  F(Ci,  C2.  C3)  det[J] 
i=l  j=l 

Cl  =AI(i) 

C2  =  AJ(j)[l-AI(i)] 

C3  =  1  -  Cl  -  C2  and 

W  =  AS(i)  H(j)  [l-AI(i)] 


The  five  sets  of  values  ( corresponding  to  the  five  point  integration  formula )  for  the 
constants  AI  (i),  AT  (j)  etc.  are  chosen  from  reference  38. 


For  the  purpose  of  numerical  computations  for  obtaining  the  stiffness  matrix  elements 
from  the  present  approach,  a  laminate  as  shown  in  Figure  5  is  considered.  Figure  5a  shows  an 
element  arrangement  for  a  quarter  plate  with  an  assumed  symmetry  about  the  x  and  y  coordinate 
axes.  Figure  5b  shows  a  set  of  more  refined  finite  element  mesh  patterns  to  be  used  later  on  in  the 
investigation.  This  is  a  typical  mesh  pattern  associated  with  triangular  finite  elements.  A  scheme 
of  nodal  numbering  is  also  shown  in  these  figures  to  obtain  an  economical  storage  capacity  in  the 
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conq>uter.  These  numbers  refer  to  the  nodes  which  have  all  the  layers  incorporated  in  than  and  thus 
the  stiffness  matrix  for  any  element  here  would  have  already  incorporated  the  stiffness  matrix 
contributions  from  all  the  layers  of  the  laminate. 


3.7  RESULTS  AND  DISCUSSION 

To  start  with,  consideration  is  given  to  laminates  with  isotropic  layer  properties.  This  is 
intended  to  check  the  validity  of  the  present  approach  by  comparison  with  known  results.  In  all  the 
following  calculations  one  global  layer  and  one  local  layer  are  considered.  From  the  mathematical 
developments  described  before,  this  gives  an  element  84  degrees  of  freedom  (the  element  containing 
both  the  local  and  global  region  as  shown  in  figure  3).  The  material  properties  chosen  are:  E  =  30  x 
10^  psi  and  G  =  11.54  x  lO^  psi.  The  geometry  chosen  is  represented  by  a  =  1"  and  total 
thickness  =  0.01". 

First,  a  square  plate  is  chosen  with  only  two  elements  in  the  quarter  plate  as  shown  in 
figure-4.  This  is  done  primarily  to  check  the  various  stages  of  calculations  in  the  procedure  such  as 
symmetry  of  stiffness  matrix,  appropriate  assembly  technique  etc.  At  the  element  level,  an 
eigenvalue  check  is  done  for  the  stiffness  matrix  which,  in  mathematical  terms,  is  given  by  the 
equation 

K  -  W  =  0 

where  K  is  the  element  stiffness  matrix,  1  is  identity  matrix,  and  A.'s  are  the  eigenvalues.  The  f  s 
(see  equation  3.42)  associated  with  X's  are  the  corresponding  eigenvectors.  The  above  solution 
must  have  at  least  six  zero  eigenvalues  corresponding  to  the  six  rigid  body  modes  for  the  structure. 
It  was  found  that  for  one  element  in  the  grid  shown  in  figure  Sa,  there  were  six  eigenvalues  very 
close  to  zero  compared  to  the  rest  of  the  eigenvalues.  It  was  also  found  that  the  diagonal  elements  of 
the  stiffness  matrix  were  nonzero,  positive  and  relatively  dominant.  A  check  was  also  done  on  the 
assembly  procedure  for  the  two  elements  shown  in  figure  5a.  Since  the  plate  is  chosen  to  be  square, 
geometric  symmetry  tests  were  also  applied.  For  example,  the  stiffness  associated  with  displacement 
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u  at  node  2  is  identical  with  that  associated  widi  displacement  v  at  node  4.  The  stiffnesses  associated 
with  w  at  nodes  2  and  4  were  found  to  be  identical.  Similar  checks  were  applied  for  symmetric 
nodes  such  as  3  and  7  and  6  and  8. 

The  next  stage  of  calculation  involved  the  actual  solution  to  the  problem  of  a  plate 
transversely  loaded  by  a  uniformly  distributed  load.  The  boundary  conditions  were  chosen  to  be 
simply  supported  on  the  edges  x  =  a  and  y  =  a.  This  implies  that  at  the  bottom  of  the  laminate  w  is 
set  to  zero  at  nodes  7, 8, 9, 6  and  3.  Also,  the  shear  stresses  '^z  ^  edges 

X  =  a  and  y  =  a  respectively.  Because  of  the  isotropic  laminate  construction  and  the  chosen 
symmetry  of  geometry,  the  following  conditions  were  also  applied. 

u  =  0  at  x  =  0,  y  =  y 

v  =  0  at  y=0,  x  =  x 

All  the  above  mentioned  conditions  were  applied  to  the  assembled  stiffness  matrix 
representing  the  entire  quarter  plate  and  the  resulting  equations  were  solved  by  a  library  subroutine. 
The  maximum  value  of  the  vertical  deflection  was  found  to  be  0.0287  qj  where  qi  is  the  intensity  of 
the  applied  load.  The  exact  value  as  given  in  standard  text  books  is 

'^max  =  0.00406  q|  a^  /  D 

where  a  is  the  side  of  the  plate  and  D  is  the  flexural  rigidity.  For  the  geometry  and  the  material 
properties  chosen,  this  deflection  is  given  as  0.02365  qi.  This  shows  an  error  of  about  21%  with 
the  finite  element  result  This  is  expected  for  the  grid  that  has  been  chosen. 

The  next  two  cases  refer  to  the  laminates  that  are  rectangular  in  planform.  For  this 
purpose,  several  grid  patterns  are  chosen  as  shown  in  figure  6.  In  the  first  case  here,  a  rectangular 
laminate  with  b/a  =  2  is  chosen.  For  this  case,  the  exact  solution  for  the  maximum  deflection  is 
given  as 
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Wmax  =  0.01013  qja^/D  =  0.00369  qj 


With  the  grids  (a)  and  (b)  in  figure  6  (i.e.,  2  and  4  elements  respectively),  the  finite 
element  results  for  were  found  to  be  0.0046  q^  and  0.00323  q^  representing  errors  of  24.7% 
and  12.4%  respectively. 

The  next  case  of  rectangular  laminate  has  an  aspect  ratio  of  4.  The  same  grid  patterns  as 
shown  in  figure  6  are  chosen  here.  However,  for  studying  the  convergence  of  results,  higher 
number  of  elements  in  the  quarter  plate  is  chosen.  These  grid  patterns  are  shown  in  figures  6c 
through  6f  in  which  the  number  of  elements  are  respectively  6, 8, 10  and  12.  The  exact  solution  for 
this  case  is 

Wmax  =  0.01282  qia^/D  =  0.000292  qj 

With  4,  6,  8,  10  and  12  elements  in  the  quarter  laminate,  the  corresponding  values  for 
maximum  deflection  were  found  to  be  0.000345  qj,  0.000325  qj,  0.000316  qj,  0.000300  qj  and 
0.000294  q^  respectively.  It  can  be  seen  that  the  results  from  the  finite  element  formulation  converge 
very  well  to  the  exact  result.  In  order  to  obtain  a  graphical  idea  about  this,  the  percentage  errors  of 
the  finite  element  results  are  plotted  against  the  number  of  elements  for  this  specific  case  with  aspect 
ratio  equal  to  4.  This  is  shown  in  figure  7.  It  can  be  seen  that  for  the  last  case  where  the  number  of 
elements  is  12,  the  error  is  about  0.7%.  These  results  establish  an  acceptable  standard  of 
convergence  rate  for  the  finite  element  stiffness  model  considered  in  this  section. 


42 


4.  GLOBAL-LOCAL  FINITE  ELEMENT  METHOD:  FRONTAL 


SOLUTION  TECHNIQUE 

In  this  section  a  set  of  problems  (bending,  uniform  tensile  and  stretch)  have  been  studied 
by  using  global-local  and  local  models.  An  inqxmant  aspect  of  this  part  is  the  use  of  frontal  method 
of  finite  element  solution.  This  is  dtme  in  older  to  reduce  the  size  of  the  structure  matrix  so  that  more 
plies  and  more  refined  mesh  can  be  studied.  Currently,  the  program  can  solve  3500  degrees  of 
freedom  or  more  on  IBM  RT  PC  without  having  any  memory  problem.  The  present  illustrative 
study  includes;  (1)  The  analysis  ofbending  effects  in  laminates  under  prescribed  uniform  transverse 
load  and  (2)  The  analysis  of  edge  effects  in  laminates  under  applied  uni-axial  uniform  tensile  stress 
or  strain.  The  basic  formulation  of  the  variational  problem  using  the  global-local  approach  has 
already  been  described  in  the  previous  section.  The  aspects  of  computer  programs  incorporating  the 
frontal  solution  technique  [42]  and  the  consequent  numerical  results  are  delineated  and  discussed  in 
this  section. 

4.1  LOCAL  DOMAIN  (MODEL  1) 

Based  on  the  reference  24,  the  simplest  assumption  consistent  with  a  realistic  stress 
analysis  for  the  in-plane  stress  component  are  assumed  to  vary  linearly  through  the  thickness  of  each 
ply.  The  substitution  of  these  stress  components  into  the  differential  equations  of  equilibrium  yields 
the  interlaminar  stress  components  in  terms  of  tractions  pj,  s^,  q  (i=l,2),  stress  and  moment 
resultants.  Finally  the  weighted  displacements  are  expressed  in  terms  of  stress  resultants  through 
constitutive  relations.  For  the  sake  of  clarity,  some  of  the  expressions  are  reproduced  here.  The 
energy  functional  is  first  written  as  [24] 

N 

%l  =  I  Z  [  aT  e  -  1/2  s  a  -  e  }  dVj^ 

V  k=l 

The  simplest  assumption  consistent  with  realistic  stress  analysis  for  the  in-plane  stress 
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component  can  be  expressed  as 


N,  12  M^z 

Ox  =  h  + 

Ny  12  M,z 

Oy  *  h  + 

N^y  12  NjyZ 

"'xy  =  h  + 

where  N^,  Ny,  N^y,  M^,  My,  M^jy  are  functions  of  x  and  y.  Obviously,  these  functions 
represent  the  force  and  moment  resultants  arising  in  plate  theory.  As  in  reference  24,  we  introduce 
the  following  stress  and  moment  resultants. 

h/2 

Nj  =  /  Oi  dz 

-h/2  i  =  1.  2,  3.  6 

h/2 

Mj  =  J  zOj  dz 

-h/2 

where  h  is  the  thickness  of  each  layer  and  N3  and  M3  are  mathematical,  not  physical,  quantities. 
We  first  consider  a  single  layer  of  thickness  h  within  the  laminate.  We  let  x  and  y  represent  the 
coordinates  in  the  midplane  of  the  layer,  which  is  bounded  by  the  planes  z=±^/2.  The  interlaminar 
stresses  o^,  x^y  at  the  top  of  the  layer  are  denoted  by  P2,  S2,  t2  respectively,  while  the 
corresponding  stresses  at  the  bottom  of  the  layer  are  designated  as  pj,  Sj,  t^. 
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In  plate  theory,  the  equations  of  equilibrium  or  balance  of  linear  momentum,  in  terms  of 
special  description  within  whole  local  region  R/  are  in  general  expressed  as 

p  “i  =  Oij.j  P  ‘‘i 

where  Oy  are  stress  components,  aj  is  the  acceleration  of  deformed  body,  bj  are  body  force  per 
unit  mass,  and  p  is  the  density  per  unit  volume.  By  neglecting  body  force  and  considering  steady 
state  only,  we  now  substitute  in-plane  stresses  along  the  values  of  the  interlaminar  stress  at  z  =  ±h/2, 
into  the  equilibrium  equations,  which  leads  to  the  following  distributions 


Oz  =  (Pl'*'  P2)  (12  z2 . 1)  +  (P2  .  Pj)  (^0^3 

4  h2  4  h 

+  3^2  (  1-  4z2  )  +  ISMj  (  2z  .  8x3  ) 

2h  h2  h2  h  h3 


Tjy  =  (*2  -  *1)  %  +  (  Sj  +$2)  (12  z2  -  1  )  +  3Vy  (  1  -  4z2  ) 

4  h2  2h  h2 


“'zx  =  ( ‘2  *  *1  )  *4  +  ( ‘l  +  ‘2>  (  12  z2  -  1  )  +  3V^  (  1  -  4z2  ) 

'1  h2  2h  h2 

where  the  shear  resultants  Vj^  and  Vy  given  by 

‘‘/2 

Vx  =  J  ^xz  *1^ 

-‘'/2 

*'/2 

Vy  =  j  Ty^  dZ 

-•'/2 
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and  are  functions  of  x  and  y  only. 


Following  theroretical  development  in  24  and  simplifing,  the  Reissner's  energy  can  be  expressed  as 


7C  —  “  7Cg 

where 

)tl  =  JJ{[V2]TlR2]T[RSlT[RND][R2][V2l+(V3)T(R3]'r[RS]T 
X  [RND]tR2][V2]  +  [V3]T(R9][V2]  +  [V3]T’(Riol  [U»] 

-  ^2  {  [V2]'^[R2l[RSl'''[RQl[RSl[R2l[V2l  +  [V3]’r[R3][RSl 
X  [RQ]T(R2][V2]  +  (V3]'r[RT][RS][R2][V2l 

+  [V2lT[R2]T'[RSl[RS][R3)T[V2]  +  [V3lT[R3]tRS] 

X  [RQ]T(R2][V3]  +  [V2]T[RT][RS][R3][V3] 

+  (V2]T(R2lT'[RS-l]T(RT]T[V3]+[V3]+T(R3]T(RS]T[RT]  x  [V2I  +  [V3]T{RP][V3]  )  dxdy 

and 

Jtg  =  J  oj  u  ds 

where 

[V2]'^  sr  [u  V  W  U*  V*  w*  w  )  =  [  U*  V*  w*  u**  ] 

[V3]T  =  (  P2  S2  S2  pi  *1  ti  ]  [VsjT  =  [  N,  Ny  N^y  My  M^^y  Vy  ] 


RND  = 


Vh 

0 

0 

0 

0 

0 

0 

0 

0 
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0 

Vh 
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RS  = 


f 

hd 

0 

0 

0 

0 

0 

0 

2dx 

0 

h_3 

0 

0 

0 

0 

0 

2dy 

0 

0 

0 

0 

0 

3 

0 

hd 

0 

0 

0 

0 

0 

2dy 

0 

0 

0 

h2d 

0 

0 

0 

R2  = 

4  dx 

0 

0 

0 

0 

d 

0 

0 

4  dy 

0 
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0 

0 

0 


0 

0 

0 


R3  = 


0  0 
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S33h  0 

10 
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S33F42 

S35F4F6  S36F4F6 

S33F3F4 

S35F4F65 

S36F4F5 

S53F4F6 

S55F6F6  S56F6F6 

S53F3F6 

S55F5F6 

S56F5F6 

S63F4F6 

S65F6F6  S66F6F6 

S63F3F6 

S65F5F6 

S66F5F6 

S33F3F4 

S35F3F6  S33F3F6 

S33F3F3 

S35F3F5 

S33F3F5 

S53F3F5 

S55F5F6  S5^5F6 

S53F3F5 

S55F5F5 

SS^SFS 

S63F3F5 

S56F5F6  S66F5F6 

S63F3F5 

S65F5F5 

S66F5F5 

Fl  = 

3/2  -  6z2/h2 

F2  = 

30z/h-  120z3/h3 

F3  = 

-V4  +  3/2Z/h+3z2/h2- 

10z3/h3 

F4  = 

*V4  +  3/2Vh+3z2/h2- 

10z3/h3 

F5  = 

-1/4+  Vh  +  3z2/h2 

F6  = 

-1/4  +  z/h  +  3z2/h2 

I3  = 

12z/h3 

I4  = 

12z/h4 

I5  = 

12z/h5 

From  above  equations,  the  funtional,  n,  has  been  defined.  Therefore  the  problem 
becomes  to  find  the  stationary  value  of  V2  and  V3  e  Hj  and  H3  respectively  such  that 

7C  =  7C/  -  Ttj  +  J  MV2dSi  +  J  NV3dS2 
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where 


H,  =  {  V2  I  V2  e  H-nCRB).  MV2  =  q,.  on  ds,  } 

H2  =  {  V3  I V3  e  H">(R3),  NV  =  q2.  on  082  } 

H'"  is  infinite  dimensional  space  dsj  and  ds2  are  the  boundary  section  in  which  V2  and 
V3  are  prescribed. 

Equilibrium  requires  that  n  be  stationary,  i.e.,  =0  where  it  must  be  recognized  that 

V2  and  V3  are  independent  variables.  Hence,  in  the  finite  element  analysis  of  on  assemblage  of 
elements  we  only  need  to  enforce  interelement  continuity  on  V2  and  V3,  which  can  readily  be 
achieved  in  the  same  way  as  in  the  isoparametric  finite  element  analysis  of  solid.  Following  the 
standard  finite  element  analysis  we  use 

q  r 

=  I  OijVi  .  V,=IXj3Vj 
i-l  j=l 

where  the  <t)i  and  Xj  are  the  interpolation  functions  and  q  and  r  are  the  number  of  nodes  of  the 
element  with  respect  to  weighted  displacement  and  interlaminar  stresses.  Using  the  fact  that  the  total 
Reissner  energy  must  be  stationary  and  taking  account  of  the  discretization,  the  final  algebraic 
equations  can  be  expressed  as 


[K] 

■  V2' 

= 

■f 

Vi 

u® 

1- 

where  the  K  is  the  stiffness  matrix  obtained  from  the  volume  integrals  of  the  energy  equations  for 
the  entire  region.  V2  and  V3  are  the  discretizing  weighted  displacements,  normal  and  shear  stress 
at  the  top  and  bottom  of  layers.  F  is  load  vector  obtained  from  the  area  integral  of  the  external 
tractions  of  whole  boundary.  U®  is  displacement  equivalent  vector  derived  from  the  interlayer 
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continuity. 


4.2  GLOBAL-LOCAL  MODEL  (MODEL  2) 

The  global-local  model  was  formulated  in  order  to  solve  the  actual  displacement  and 
interlaminar  stresses  at  the  same  time.  The  energy  equation  for  the  global  regions  and  local  layer  can 
be  expressed  as 

7C  =  J  WdV  +  J  (1/2  Oij  (Ujj  +Uy)  -  W]dV  -  J  TiUjds 
V  ''  s 

where  W  =  W  (Uj,  Cy).  W  =  W  (oy,  Cy)  and  W  and  W  are  strain  energy  density  functions,  the  first 
in  terms  of  displacements  Uj  and  ejj.  the  expansional  strain  components,  and  the  second  in  terms  of 
stresses  Ojj  and  ejj. 

The  energy  equaticm  can  be  written  as 

K  =  llg  +  71/  -  Ks 

where 

TCg  =  J  [U]  {  q'*'  A  qi  +  qTB  qj  +  qT  B  q2  +  qT’E  q2 

Ag  +  qTEqj  +  qTFj  q3  +  q’’’ Ej  qj  +  qTpi  qj 

+  qTHiq3  }(V)dxdy 
N 

7C/  =  /  Z  {  [U]'''’(D6]T[rs][RQJ(RSJ[D5][V]  +  (UjjTjRjjTjRslTjRQjjDgjjV] 

A/  k=l  +  [V]T(D6][RsF(RQJIRS]  [R3][U3]  +  [U3](R3](RS]RQ](RS](R3]  (V3] 

+  (U3]T(RTl(RS)[R3l(V3]  +  [V3]T(R3)T[rsiTirt][U31 

+  IUl'^lD6]T[RSjTiRT]IV3]  +  (VjKRTKRSJIDgKU] 

+  [U3]T[RP][V3]  dxdy 

Kg  =  /  fg  u  d  Sg  +  /  fj  u  d  Sj 

Sg  S/ 
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where  the  and  A/  represent  the  area  of  the  global  and  local  region  in  the  x-y  plane  respectively. 
The  Sjj  and  S/  indicate  the  external  boundaries  of  the  global  and  local  regions  respectively.  The  f^ 
and  f]  stand  for  the  external  forces  acting  at  the  global  and  local  regions  and 
Qi  =  Dj  D4  D5-1 
02  =  D2  D4 
03  =  D3  D4  D5-1 
Dfi  =  R2  D5-1 


52 


54 


(  A,  B,  E  )  = 


h/2 

J  (  1,  z.  z2  )  C  dz 
-h/2 

(Ei.Fi)  =  J  (J)(z2.  z3)c3  dz 

-h/2 

h/2 

Hi  =  J  (.5)(z4)M22  dz 

-h/2 


and  [R2I.  [R3].  [RS],  [RQ],  [RP],  and  [RT]  are  defined  as  local  model  in  the  previous  section. 


From  above  equations,  the  functional  Jiof  71^.71/  and  tc,  are  defined.  Therefore  the 
problem  becomes  to  find  the  stationary  value  of  U,  U3  e  and  H2  respectively  such  that 


where 


7t  =  Ttg  +  TCj  -  TCj  + 


J  MVdSj  + 


dSj 


J  NV3  dS2 
dS2 


Hi  =  {  VIVe  (R3).MV  =  qi  on  dSi  } 

H2  =  {  V3  I  V3  e  H™  (r3).  MV3  =  qj  on  aS2  } 


H*"  is  infinite  dimensional  space,  dSj  and  as2  are  the  boundary  section  in  which  U  and 
U3  are  prescribed. 


Equilibrium  requires  the  7C  is  stationary,  i.e.  dTt  =  0  where  it  must  be  recognized  that  U 
and  U3  are  independent  variables.  Hence  in  the  finite  element  analysis  of  an  assemblage  of  elements 
we  only  need  to  enforce  interelement  continuity  on  U  and  U3,  which  can  readily  be  achieved  in  the 
same  way  as  in  isoparametric  finite  element  analysis  of  solid.  Following  by  standard  finite  element 
analysis  we  choose 
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q  r 

U  =  Z  <l»i  Ui  ,  U3=  I  Xj  sUj 

i=l  j=l 

where  the  *14  and  Xj  are  the  interpolation  functions,  and  q  and  r  are  the  number  nodes  of  the 
element  with  respect  to  weighted  displacement  and  interlaminar  stresses. 

Using  the  fact  that  the  total  energy  must  be  stationary  and  taking  account  of  the 
discretization,  the  final  algebraic  equation  can  be  expressed  as 

=  [F] 

where  the  K  is  the  stiffness  matrix  obtained  from  the  volume  integrals  of  the  energy  equations  for 
the  entire  region.  U,  U3  are  the  discretizing  actual  displacements,  normal  and  shear  stress  at  the  top 
and  bottom  of  layers.  F  is  load  vecttn*  obtained  from  the  area  integral  of  the  external  forces  of  whole 
boundary. 

The  computer  code  "GLFEM",  based  on  the  mixed  variational  finite  element  formulation 
with  frontal  solver  have  been  successfully  developed.  A  transverse  loading  (bending)  of  simple 
supported  case  have  been  studied  to  test  the  numerical  solution  and  convergence.  Also  a 
unidirectional  tension  load  applied  at  the  edge  of  the  plate  have  been  studied.  The  flow  chart  of  the 
program  is  given  in  Appendix  A,  and  in  a  brief  description  of  the  program,  "GLFEM",  is  as  follows: 

i)  Data  files  have  been  opened  by  the  Index  of  NFIA,  NFIB,  NFDA,  NFOB,  NFOC, 
NFOD  and  MY,  the  first  two  Index  are  input  files  and  others  are  output  files.  If  number 
of  variable  is  more  thatn  3,000,  then  other  output  file  is  recommended  to  open.  All 
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the  data  files  are  open  at  beginning  of  main  program  and  close  at  end  of  main  program, 
for  more  detail  see  appendix  A. 


ii)  All  the  basic  data  is  input  or  generated  in  the  INPUT  subroutine.  The  data  include 
number  of  total  element,  node  coordinates,  node  numbers,  element  connectivity,  control 
index,  material  properties  for  each  layer,  integrating  data,  boundary  condition  data  and 
control  data  of  pre-fiontal  part. 

iii)  The  element  data  based  on  two  dimension  model  can  be  directly  read  from  data  files 
generated  from  pre-processor  program.  By  using  convert  subroutine  those  data  are 
converted  into  3  dimensional  data.  The  basic  FRONTAL  routine  data  are  generated  in 
INPUT  subroutine. 

iv)  The  smeared  (or  effective  moduli)  for  both  global  and  local  region  are  calculated  by 
GLOBAQ  and  LOCAQ  subroutine  respectively.  The  pre-data  of  global  and  local 
regions  are  obtained  by  SMAT,  USBAR  and  ALLOT  subroutine. 

v)  The  shape  functions  and  their  derivatives  for  both  linear  and  quadratic  functions  in 
triangular  and  rectangular  elements  are  set  up  in  SFR2,  SHAPW,  SHAPQ2  and 
SHAPL  subroutine.  The  matrices  of  the  transformation  are  evaluated  in  JACOB2 
subroutine.  The  Gauss-Radau  data  is  input  in  GAUSPT  subroutine. 


vi)  The  stiffness  matrix  of  global  region  is  calculated  in  ESTIFPS  subroutine.  The 
stiffness  matrix  of  local  region  is  evaluated  from  DBLOCA  subroutine.  The  predata 
values  are  calculated  by  PREPBM,  DISWET  subroutine,  the  inverse  of  the  matrix  is 
executed  by  SOLVE  subroutine. 

vii)  The  load  vector  can  be  composed  by  concentrated  and  body  forces,  edge  distriution 
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forces  and  distributed  bending  forces,  which  are  coded  by  LOADPS,  BODY, 
EDGEF  and  LOADPB  subroutine  respectively. 

viii)  The  assembly  of  the  stiffness  matrix  and  elimination  of  the  variables  for  the  global 
and  local  domains  arc  done  by  FRONTAL  subroutine.  The  standard  frontal 
procedure  are  followed  (for  more  detail,  see  reference  42). 

4.3  EXAMPLE  PROBLEMS  AND  NUMERICAL  RESULTS 

The  theoretical  formulation  is  made  for  the  analysis  of  composite  structures  with  any 
geometric  shape  and  size.  Three  cases  of  boundary  value  problems  have  been  investigated.  Case  a 
deals  with  the  bending  of  isotropic  thick  and  thin  rectangular  plates  with  simply  supported  edges  due 
to  a  uniformly  distributed  transverse  load.  Case  b  deals  with  the  stress  analysis  of  composite 
laminate  under  the  influence  of  an  applied  in-plane  unidirectional  stress  at  the  edges  (x=constant)  of 
the  laminate.  Case  c  deals  with  the  stress  analysis  of  a  composite  laminate  under  the  influence  of  an 
applied  in-plane  strain  at  x=constant.  The  results  are  presented  in  graphical  form.  A  comparison  is 
made  with  existing  results  and  show  a  good  agreement. 

Case  a:  Bendinecase 

For  tesing  the  computer  program  written  for  the  model  2,  the  solution  to  the  problem  of  a 
plate  transversely  loaded  by  a  uniformly  distributed  load  with  simply  supported  edges  was  obtained. 
The  laminate  geometry  is  shown  in  figure  8.  Because  of  the  isotropic  laminate  construction  and  the 
chosen  symmetric  geometry  (see  figure  6a),  the  following  conditions  were  applied: 

u=0atx  =  0  v  =  0aty  =  0 

The  boundary  conditions  were  chosen  to  be  simply  supported  on  the  edges  where  x  =  a,  y  =  b. 
This  implies  that  at  the  edge  of  the  laminate  the  transverse  deflection  W  is  set  to  zero  for  all  nodes. 
Also,  the  shear  stress  components  and  x^y  are  set  to  zero  on  edge  x  =  a  and  y  =  b 
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respecti  /ely.  All  the  above  mentioned  conditions  were  applied  to  the  assembled  stiffness  matrix 
representing  the  actual  structure  plate.  The  above  solution  of  square  plate  as  given  in  standard  text 
books  is 


Wmax  =  0.00406  qa4/D 

where  q  is  the  intensity  of  the  applied  loads,  a  is  the  side  of  the  plate  and  D  is  the  flexural 
rigidity.  Material  properties  in  the  plane  of  elastic  symmetry  of  each  layer  are  given  by 

El  “  Et  —  E^  ~  30  X  lO^psi,  Glx  ~  Gf y  =  ®tz  ~  1.154  x  lO^psi 
Vlt  =  ^LZ  -  ^TZ  -  0.3, 

where  L,  T  and  Z  refer  to  fiber,  transverse,  and  thickness  directions  respectively.  The  exact 
maximum  deflection  for  square  plate  is  given  as  0.02365  q  and  0.002956  q  for  plate  thicknesses 
equal  to  .01  and  .02  inch,  respectively.  The  exact  maximum  deflection  for  rectangular  plate  with  a/b 
=  2  is  given  as  0.00369  q  and  0.000461  q  for  the  plate  thicknesses  equal  to  0.01  and  0.02  inch 
respectively. 

The  triangular  elements  and  nodes  are  shown  in  figure  3.  Figure  5  shows  the  grid 
patterns  used  in  the  model  (for  a  square  plan).  With  grid  patterns  1,2,  3  and  4  in  the  quarter  laminate, 
the  corresponding  values  for  maximum  deflection  listed  in  these  tables.  It  can  be  seen  that  the 
results  from  the  global  local  finite  element  model  converge  very  well  to  the  exact  value.  In  order  to 
obtain  a  graphical  idea  about  this  aspect,  the  displacement  parameter  W  calculated  by  the  finite 
element  method  and  given  in  Tables  4.1  and  4.2  are  plotted  in  figures  9  and  10  against  the  number 
of  grids  of  the  square  and  rectangular  plates  for  both  plate  thicknesses. 
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1a  Thickness  =  .01  INCH 


grid  number 

W  FEM  RESULTS 

EXACT  SOLUTION 

1 

.0185  q 

.02365  q 

2 

.0189  q 

3 

.0196  q 

4 

.0206  q 

1b.  Thickness  =  .02  INCH 


W  FEM  RESULTS 
.00232  q 
.00240  q 
.00258  q 
.00274  q 


Table  4.1  Rectangular  Plate  Solution 


2a  Thickness  =  .01  INCH 


grid  number 

W  FEM  RESULTS 

EXACT  SOLUTION 

1 

.00297  q 

.00369  q 

2 

.00306  q 

3 

.00324  q 

4 

.00342  q 

2b.  Thickness  =  .02  INCH 


grid  number 

W  FEM  RESULTS 

EXACT  SOLUTION 

1 

.000373  q 

.0004609  q 

2 

.000396  q 

3 

.000427  q 

4 

.000444  q 

Table  4.2  Square  Plate  Solution 
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By  comparing  the  figures  9  and  10  in  numerical  sense,  the  error  percentage  change  from 
7.4%  to  3.6%  as  the  total  plate  thickness  changes  from  .01  inch  to  .02  inch  for  the  rectangular  plate 
using  grid  4.  Similarly  for  square  plate,  the  error  percentage  changes  from  12.6%  to  7.26%.  The 
results  in  figures  9  and  10  show  that  the  converging  speed  and  error  percentage  of  rectangular  plate 
are  better  than  those  of  square  plate. 

Case  b:  Aonlied  Tensile  Stress 

This  section  shows  the  solution  for  the  problem  of  a  finite  plate  subjected  to  an  in-plane 
tensile  load  Cx  =  1  at  the  edge,  x  =  a  by  using  models  1  and  2.  The  laminate  geometry  is  shown  in 
figure  8.  Two  cases  of  geometric  conrigurations,  the  ratio  of  a/b,  are  studied  using  the  first  and 
second  model.  Also,  the  solutions  for  the  thin  and  thick  finite  plate  are  compared.  A 
symmetric  laminate  with  the  stacking  of  layers  as  0®,  90®,  90®,  0®  will  be  examined  in  this  study. 
Comprehensive  results  based  on  the  reference  24  will  be  employed  to  compare  specific  results  given 
by  the  present  theory.  The  layers  are  of  equal  thickness  h,  the  laminate  width  is  2b  =  80h  or  2b  = 
16h  and  material  properties  in  the  plane  of  elastic  symmetry  of  each  layer  are  given  by 

Ell  =  20  X  10^ psi ;  E22  =  E33  =  2. 1  x  10^ psi 

Gi2  =  Gi3  =  G23  =  .85  X  10^  psi  V12  =  i>13  =  V23  =  21 

where  1,2  and  3  refer  to  the  fiber,  transverse,  and  thickness  directions  respectively,  and  t)i2  for 
example,  is  the  Poisson  ratio  measuring  strain  in  the  transverse  direction  due  to  uniaxial  tension  in 
the  fiber  direction. 

Because  of  the  chosen  symmetry  of  geometry  in  Figure  5a,  the  following  conditions  were 

applied 

For  model  1:  u,  u*  =  0  at  x  =  0  v,  v*  =  o  at  y  =  0 

For  model  2:  u‘,  u**  =  0  at  x  =  0  v‘,  v*’  =  o  at  y  =  0 
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where  parameters  in  model  1  are  weighted  displacements  in  x  and  y  direction  and  are  defined  in 
references  24, 25  and  the  superscripts  t  and  b  denote  the  respective  variable  at  top  and  bottom  of  the 
layer,  respectively.  For  applied  in-plane  tensile  load  case,  the  boundary  conditions  of  the  finite  plate 
for  both  formulations  are  as  follows: 

The  shear  stress  components  X2x  'fzy  are  set  to  zero  at  x  =  a  and  y  =  b  respectively. 
At  mid-surface,  the  transverse  deflection  W  and  shear  stress  components  x^y  and  Xj^  are  set  to  7.ero 
everywhere.  All  the  above  mentioned  conditions  were  applied  to  the  assembled  stiffness  matrix 
representing  the  actual  structure  plate. 

The  triangular  elements  and  nodes  are  shown  in  Fig.  3  and  different  grid  patterns  are 
shown  in  figure  5b.  Figures  11  to  16  show  the  convergence  study  of  both  interlaminar  stress 
components  and  weighted  displacement  for  thick  plates  of  aspect  ratio  a/b=5.  The  applied  load  is 
constant  in-plane  tensile  stress  at  the  edge  of  x=constant.  The  results  obtained  from  model  1  have  the 
same  trend  compared  with  reference  24,  but  magnitudes  are  slightly  different.  That  is  because  the 
applied  boundary  condition  in  the  present  formulation  is  Ox  =  1  at  x  =  a  while  in  reference  24  is 
6^=  e.  The  results  of  the  stresses  from  the  model  2  are  almost  equal  to  zero  at  free  edge  even  though 
the  actual  displacements  are  reasonable.  Figure  17  shows  the  comparing  interlaminar  stress  of 
different  geometric  parameter  for  applied  tensile  load. 

Case  c:  Uniform  Applied  Strain 

In  this  section  the  response  of  the  laminate  under  the  influence  of  the  unidirectional 
constant  strain  applied  at  the  edge  of  the  plate  has  been  studied.  This  case  is  simulated  by  applying 
the  equivalent  biaxial  stress  obtained  by  lamination  theory.  The  convergence  study  of  both 
interlaminar  stress  and  weighted  displacement  of  the  thin  and  thick  plate  is  similar  to  the  one  with 
uniform  tensile  load  applied  at  the  edge.  Figures  18  to  22  show  the  convergence  pattern  of 
interlaminar  stresses  and  weighted  displacement  for  thick/thin  plate  of  aspect  ratio  a/b=5.  The 
results  of  interlaminar  stresses  obtained  from  model  1  for  constant  stain  case,  have  the  same  trend 
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but  different  magnitude  as  reference  24.  Figures  23  to  25  show  the  comparison  of  interlaminar 
stress  distribution  between  the  existing  and  present  results.  The  differences  in  magnitude  between 
the  current  iiKxlel  and  those  of  reference  24  may  be  explained  by  the  following  points: 


1 .  The  present  is  3D  model  and  solution  is  obtained  fcx*  finite  length  plate. 

2.  The  meshes  are  not  fine  enough  near  the  edge  and/or  more  sublayers  should  be  used 
in  order  to  get  better  accuracy  in  results. 

Figures  23  to  25  also  show  the  comparison  of  interlaminar  stress  of  thick  (h=.25  in.)  and  thin 
(h=.05  in.)  plate  for  both  strain  load  with  respect  to  aspect  ratio  a/b=5. 

4.4  SUMMARY  AND  CONCLUDING  REMARKS 

The  convergence  of  current  model  is  good  as  shown  in  the  pictures.  The  results  of 
interlaminar  stress  components  and  weighted  displacement  for  both  constant  tensile  and  strain  load 
from  the  Model  1  have  the  same  trend  but  different  magnitude  compared  with  those  of  reference  24. 
The  differences  in  magnitude  between  the  current  model  and  existing  results  may  be  explained  as 
follows: 


1.  The  present  is  3D  model  and  solution  is  obtained  for  finite  length  plate. 


2.  The  meshes  are  not  fine  enough  near  the  edge  or  more  sublayers  should  be  used 
in  order  to  get  accurate  results  especially  for  the  thin  plate. 


3.  Applied  boundary  conditions  are  different. 


4.  Geometric  parameters  are  different. 


The  results  of  interlaminar  stresses  obtained  from  Model  2  are  almost  zero,  even  though  the 
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displacement  trend  makes  sense.  The  reason  may  be  explained  as  follows: 


1.  The  meshes  are  not  fine  enough  near  the  edge.  For  more  accurate  results,  more 
refined  mesh  and/or  more  sub-layers  should  be  used. 

2.  The  transformation  used  to  convert  the  weighted  displacements  into  actual 
displacements  may  introduce  additional  applied  loads,  even  though  the  displacement 
makes  sense. 

3.  The  global  region  may  be  needed  to  include  the  interlaminar  stresses  as  well  as 
displacements. 

For  the  thin  plate,  the  results  in  Hgures  23  to  25  show  that  the  shear  stresses  are  more 
significant  than  the  normal  stress  as  compared  to  thick  plate  solution.  It  should  be  pointed  out  that 
because  the  mesh  near  the  edge  is  not  very  fine,  the  models  do  not  give  prominent  peak  stresses  in 
the  presented  thin  plates  solutions.  This  is  because  of  the  limited  current  hard  disk  space.  A  more 
refined  mesh  will  be  studied  in  future  investigations.  Convergence  study,  as  shown  in  pictures, 
indicates  that  the  results  for  constant  strain  applied  load  case  are  more  stable  than  those  for  the 
contant  tensile  load  case.  The  results  by  using  the  constant  tensile  load  show  slight  oscillation. 
These  oscillations  reduce  very  significantly  as  the  clement  number  increases  near  the  edge.  The 
aspect  ratio  (a/b)  may  also  have  effect  on  the  oscillation  for  the  tensile  load  case. 

The  current  code  uses  frontal  solver  and  helps  reduce  memory  space  requirements.  For 
example  the  program  can  handle  3500  degrees  of  freedom  on  IBM  RT  PC  only  using  about  400 
maximum  frontal  size,  but  the  execution  time  is  increased.  The  benefit  of  using  frontal  solution  will 
increase  with  the  solution  of  laminates  with  large  number  of  layers.  The  future  studies  include, 
modification  of  the  computational  procedure  and  the  application  of  this  model  for  investigating 
laminates  with  flaws  and  through  the  thickness  holes.  Also  the  curing  stress  study  will  be  included. 
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5.  STIFFNESS  AND  STRENGTH  DETERMINATION  OF  MULTI- 


ORIENTED,  COATED  FIBER  COMPOSITES 

In  this  part  of  the  report,  a  coated  continuous  fiber  composite  has  been  analyzed  using  a 
three-phase  concentric  cylinder  model.  Effective  thermoelastic  properties  have  been  determined  and 
the  solution  for  the  stress  distribution  derived  under  a  uniform  three  dimensional  mechanical  and/or 
hygrothermal  loading.  First,  the  theoretical  model  is  described  and  the  numerical  results  are 
compared  with  some  of  the  closed  form  solutions  already  available  in  the  literature.  Next,  details  of 
the  computer  program,  called  NDSANDS,  which  was  developed  to  implement  the  numerical 
alogorithm,  are  discussed.  Finally,some  numerical  results  are  presented  for  the  effective 
thermoelastic  moduli  and  the  micro  stress  distribution  under  a  uniform  temperamre  change  for  a 
three-dimensional  coated  fibrous  composite. 

5.1  THE  COMPOSITE  MODEL 

A  coated ,  continuous  fiber  reinforced  composite  is  modeled  by  a  representative  volume 
element  composed  of  N  concentric ,  circular  cylinder  elements,  in  which  the  innermost  cylinder  is  the 
fiber,  the  next  ring  is  the  coating  and  the  outer  ring  is  the  matrix  as  shown  in  figure  26.  Let  us  further 
denote  the  composite  volume  in  between  the  elements  as  the  interstitial  matrix  region.  Both  the  matrix 
in  the  composite  cylinder  as  well  as  in  the  interstitial  region,  in  turn,  could  be  reinforced  by  particles. 

0)  a) 

Fitch  element  orientation ,  j ,  is  defined  via  the  two  cylindrical  angles  G  and 

with  respect  to  a  fixed  Xj  -  X2  -  X3  coordinate  system.  The  local  element  cartesian  coordinate 
system  is  represented  by  X]  -  X2  -  X3  .  It  should  be  further  noted  that  the  local  fiber  axis,  X| , 
coincides  with  the  z  -  coordinate  in  the  local  cylindrical  system,  whereas,  X2  -  X3  is  the 
transverse  r  -  0  plane. 

0)  (j) 

If  we  denote  Qu  =  cos(  Xj^ ,  X] ),  then  from  geometry  in  figure  26,  we  have 
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(j) 

Old  = 


(5.1) 

5.2  MODEL  ASSUMPTIONS  AND  METHOD  OF  ANALYSIS 

The  general  definition  of  the  effective  elastic  moduli  of  heterogeneous  materials  has  been 
discussed  in  [43]  and  [44].  For  the  purpose  of  self-consistency,  a  short  discussion,  specific  to  the 
problem  here  treated,  will  now  be  given. 

In  nber-reinforced  materials,  the  ratio  of  length  to  fiber  diameter  is  usually  very  large. 
Accordingly,  fiber  end  conditions  will  only  be  considered  here  in  the  Saint  Venant  sense.  The  fiber, 
coating  and  the  matrix  are  assumed  to  be  linearly  elastic,  homogeneous  and  perfectly  bonded.  In 
general,  the  constituent  materials  may  have  transversely  isotropic  elastic  and  thermal  expansion 
coefficients. 


(j) 

cosf2 

0 

-  sin  Q 
0 


0)  (j) 

sinQ  cos(|> 

(j)  0) 

cosQ  cos<|> 

O') 

-  sin<{> 


0) 

sin  Q 

0) 

cos  Q 

0) 

COS(> 


(j) 

sm<t> 

0) 
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Let  the  composite  material  volume  be  subjected  to  a  set  of  boundary  conditions  of  the 


form 


Ui(s)  *  Cij  Xj 


or  Ti(s)  =  Ojjnj 


(5.2a,5.2b) 


where  nj  is  the  unit  outward  normal  vector  on  the  boundary  surface  S,  xj  are  the  cartesian 

0  0 

coordinates  of  that  surface  and  e^j  and  Ojj  are  constants. 

For  (5.2)  prescribed,  it  can  be  shown  that 


£,]=  ^j=  const. 


or 


Oii  = 
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0 

Oij  =  const 


(5.3a,5.3b) 


respectively,  where  an  overbar  denotes  the  average  value  over  the  whole  volume.  When 

0 

displacements  are  presciibed,  the  average  strains  are  ejj  and  Cy  have  to  be  found  and  for 

0 

prescribed  traction,  the  average  stresses  are  Cyand  Cij  need  to  be  determined. 


The  specimen  is  further  assumed  to  be  noacroscopically  homogeneous,  by  which  is  meant 
that  for  (5.2)  prescribed,  strain  and  stress  averages  taken  over  large  enough  subregions  of  the 
specimen  are  the  same  for  any  such  subregion.  Such  a  subregion  will  be  referred  to  as  a 
representative  volume  element  (RVE). 

To  facilitate  the  analysis  we  next  introduce  an  equivalent  homogeneous  medium  (figure 
27),  having  the  effective  composite  properties,  as  a  comparison  material  and  assume  that 
the 

a) 

displacements  or  traction  acting  on  the  boundary  of  each  composite  cylinder  element ,  Sc ,  can  be 
approximated  by  those  on  the  comparison  material ,  i.e.. 


Ci)  0  0)  (j) 

Uj  =  Eik  Xjc  on  Sc 


(j)  0  (j)  (j) 

or  Tj  =  Ojjc  njc  on  Sp 


(5.4a.5.4b) 


according  to  whether  (5.2a)  or  (5.2b)  is  prescribed  and  no  summation  on  the  superscripts  is 
implied.  Further,  equation  (5.4b)  applies  on  the  radial  boundary  of  the  composite  cylinder 
elements  while  on  the  ends  the  following  resultants  need  to  be  specified: 


2n  r3 


r  r 

=  OzzT  dr 

J  « 


d0  ; 


0  0 


2  n  r3 


dr  d0 


0  0 


2  It  To 

/•  /•  *5 


•  ^ 
J 

0  0 


r^  sin  0  dr  d0 


2  It  r3 


M3  = 


Ozz  cos  0  dr  d0 


0  0 
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where  Fj,  is  the  axial  force,  T  the  torque,  and  M2  and  M3  the  resultant  moinrat  about  the  X2  and  X3 
axis,  respectively.  The  resultant  forces  in  the  X2  and  X3  directions,  F2  and  F3  respectively,  are 
non-zero  but  are  not  independent  boundary  conditions.  Thus,  under  boundary  conditions  leading  to  a 
constant  strain  or  stress  field  within  a  homogeneous  body,  the  stress  and  the  displacement  inside  the 
composite  cylinder  element,  can  now  be  evaluated  as  described  in  the  following  section. 

5.3  DETERMINATION  OF  STRESS  AND  DISPLACEMENT  FIELDS 

If  the  composite  volume  is  now  subjected  to  boundary  conditions  given  by  (5.2a)  or 
(5.2b),  then,  within  each  of  the  N  elements,  we  have  three  displacement  fields  (in  materials  (1),  (2) 
and  (3)  or  fiber,  coating  and  matrix,  respectively ).  The  form  of  the  governing  field  equations  and 
boundary  conditions  lead  to  a  general  solution  of  the  type. 

(j.p)  2  (j.p)  fcosne 

ui  =  Z  Fin(r,  z)  \  sin  n0  (5.6) 

n  =  0 

where  j  =  1.  2 . N  ;  p  =  1, 2,  3  . 

Using  the  strain  displacement  equations  and  the  stress-strain  relations  for  a  transversely 
isotropic  constituent,  the  stress  field  is  expressed  as 

(j.  p)  2  (j.  p)  f  cos  n0 

Oi  =  Z  Hin(r,  z)  l^sin  n0  (5.7) 

n  =  0 

(j.  P)  (j.  P) 

where  the  general  form  of  the  functions  Fin(r,  z)  and  Hin(r,  z)  in  equations  (5.6)  and 
(5.7),  respectively,  are  given  in  Ref  [45]  and  the  specific  constants  are  to  be  evaluated  by  the 

following  boundary  /  interface  conditions: 

i)  Displacements  (eq  5.4a)  or  traction  (eq  5.4b)  and  end  resultants  (eqs  5.5)  are  prescribed 
according  to  (5.2a)  or  (5.2b),  respectively,  at  the  boundary  of  the  composite  cylinder 
assemblage  ;  Equilibrium  of  the  entire  body,  however,  imposes  certain  implied 
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connections  among  the  traction  ccnnponents  [46]; 

ii)  Traction  and  displacements  must  be  continuous  across  the  interfaces;  and 

iii)  Stresses  and  displacements  must  be  bounded  at  the  origin,  r  =  0. 

For  convenience,  the  terms  involving  rigid  body  motion  can  be  further  neglected.  The 
expressions  for  the  constants  are  available  in  detail  in  Refs  [46, 47]. 

5.4  THE  COMPOSITE  RESPONSE 

The  cotrqrosite  stress,  or  strain  Ey,  can  now  be  determined  by  volume 
averaging  the  stress  or  strain  field  over  the  oxistituents,  namely,  the  composite  cylinder  elements  and 

the  interstitial  matrix,  respectively.  The  stress-strain  relation  for  the  composite  now  takes  the  form 

®kl  “  C|dinn  (  £mn  "  ®mn)  ^  “  Sumn  ®mn  +  ®kl  >  (5.8a,5.8b) 

where  Cumn  is  the  effective  stiffness,  Su^n  is  the  effective  compliance  and  e^i  is  the 
effective  expansional  (non-mechanical)  strain  of  the  composite. 

To  evaluate  the  effective  elastic  moduli,  we  set  the  expansional  strain  components 
identically  equal  to  zero,  i.e., 

(i.  p) 

Cy  =  0  (for  all  j.  p )  (5.9) 

With  (5.9),  the  stress-strain  relation  for  the  composite,  eqs.  (5.8),  therefore  reduce  to 
(using  contracted  notation) 

—  C|f(  £|  or  £jj  =  Sjfi  0[  (  k,  I  =  1,  2,...  6)  (5. 10a,5. 10b) 
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By  setting  each  strain  (or  stress ) component  equal  to  one  individually,  while  all 


others  are  zero,  we  will  respectively  obtain  the  1*  column  of  the  Cu  (  or  Su)  matrix. 
The  composite  engineering  constants  can  now  be  defined  in  terms  of  the  elastic  stiffnesses  (or 

compliances).  Finally,  the  effective  expansional  strains  can  be  determined  by  subsituting  in  eqs. 
(5.8). 


5.5  INITIAL  FAILURE  CRITERION 

The  problem  of  the  analysis  of  failure  of  C(»nposite  materials  is  by  an  order  of  magnitude 
more  difficult  than  the  problem  of  physical  property  prediction  which  has  been  discussed  until  now. 
When  a  composite  specimen  is  subjected  to  increasing  load  and/or  temperature,  microfailure  will 
develop  at  some  stage.  These  may  be  in  the  form  of  matrix  cracks,  fiber  ruptures,  interface 
separation  and  local  plastification.  As  loading  continues,  they  will  multiply  and  ultimately  merge  to 
produce  catastrophic  failure. 

The  criteria  for  strength  or  failure  are  often  expressed  in  terms  of  stresses,  but  such 
criteria  do  not  necessarily  refer  only  to  a  state  of  complete  rupture  of  the  materials.  Failure  criteria 
may,  in  fact,  refer  to  the  initial  events,  as  evidenced  in  yielding,  that  are  the  precursors  of  ultimate 
failure. 


One  such  simple  criteria,  which  is  widely  used  to  predict  the  initial  or  Erst  failure  of 
composite  materials,  is  the  so  called  maximum  stress  criteria.  This  criteria  is  applied  by  calculating 
the  strength/stress  ratio  for  each  stress  component  The  sign  of  the  normal  stress  determines  if 
tensile  or  compressive  strength  should  be  used  (primed  quantities  mean  compressive).  The  lowest 
strength  ratio  among  the  following  three  equations  determines  the  ratio  that  controls  the  failure  of  the 
composite. 


Rx  =  X  Ox>0,  or 


R*'  =  _2C_  if  a^<0 

10,1 


(5.11) 

contd. 
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or 


(5.11) 


Ry  =  Y  if  Oy  >  0, 


R,  =_S_ 
lasi 


R'  =  _Y1  ifO„<0 
I  Gy  I 


where  X  =  Longitudinal  tensile  strength 

X'=  Longitudinal  compressive  strength 
Y  =  Transverse  tensile  strength 
Y'  =  Transverse  compressive  strength 
S  =  Longitudinal  shear  strength 


When  the  applied  stress  component  is  unity,  the  resulting  strength  ratio  is  the  strength. 

Most  of  the  other  failure  theories  of  isotropic  materials  are  applied  on  the  principal 
stresses  ( Oi,  02  and  03 )  within  the  structure.  The  state  of  stress  at  a  point  is  first  transformed  (via 
Mohr's  circle  for  instance)  to  obtain  the  principal  stresses. 

For  ductile  materials,  Tresca  proposed  the  following  yield  criteria 

max  {  0.5  I  Oi  -  02 1 , 0.5  I O2  -  O31 , 0.5  I  03  -  Ojl  )  =  Xy  (5.12) 

where  Xy  is  the  yield  limit  in  simple  shear,  which  is,  for  ductile  materials,  set  equal  to  the  shearing 
stress  at  yield  in  simple  tension  or  compression. 

According  to  the  maximum  normal  stress  theory,  a  brittle  material  fails  when  any  of  the 
principal  stresses  reaches  the  ultimate  value,  i.e., 

max  {  I  Oi  I ,  I  02  I ,  I  03  I  }  =  Ou  (5.13) 
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The  maximum  stress  criteria  is  often  used  because  of  its  simplicity,  but  it  has 
significance  for  one-dimensional  state,  and  its  use  in  multidimensional  stress  state  must  be  made  with 
caution.  Again  we  should  avoid  extending  die  sin^listic  failure  modes  based  on  maximum  stress 
components  to  fiber,  matrix  and  inteifacial  failure  modes  because  the  micromechanics  of  failure  is 
highly  coupled  and  does  not  reduce  to  clearly  defined  modes  of  two  or  three  types. 

5.6  VERinCATION  OF  THE  THEORETICAL  MODEL 

To  p’ace  the  present  theory  in  proper  perspective  and  to  test  the  validity  of  the  analysis, 
we  numerically  simulated  some  trial  cases  and  made  comparisons  with  some  of  the  already  published 
exact  solutions.  The  reported  results  are  based  upon  displacement  boundary  formulation  since  this 
has  traditionally  led  to  better  agreement  with  experimental  measurements  of  composite  moduli. 

a.  Effective  ModuU 

Expressions  and  bounds  for  the  five  effective  elastic  moduli  of  a  unidirectional  fiber 
composite,  consisting  of  transversely  isotropic  fibers  and  matrix,  had  been  derived  by  Hashin  [48] 
on  the  basis  of  analogies  between  isotropic  and  transversely  isotropic  elasticity  equations.  Using  the 
properties  of  a  graphite/epoxy  system,  the  bounds,  along  with  the  calculated  results,  at  various  fiber 
volume  fraction,  c  ,  are  listed  in  Table  5.1  below.  ( values  within  paranthesis  indicate  lower  and 
upper  bound  respectively) 
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Table  5.1:  Comparison  of  elastic  moduli 

(Hashin's  results  and  present  analysis) 


c  EAXlO^MPa 

Va 

G^xlO^MPa 

Kx  103  MPa 

GjX  103  MPa 

58.116 

.32407 

1.3804 

4.5697 

1.4931 

0.16  (58.116, 

(.3217, 

(1.3804, 

(4.5684, 

(1.4766, 

58.131) 

.32407) 

1.3855) 

4.5972) 

1.5141) 

88.8597 

.3098 

1.4404 

4.7619 

1.6407 

0.25  (88.8597, 

(.3067, 

(1.4404, 

(4.7530, 

(1.6045, 

88.8807) 

.3098) 

1.4478) 

4.7955) 

1.6622) 

126.433 

.2927 

1.5175 

5.008 

1.8465 

0.36  (126.433, 

(.2891, 

(1.5175, 

(4.995, 

(1.7817, 

126.458) 

.2927) 

1.5271) 

5.051) 

1.8622) 

Four  out  of  the  five  effective  elastic  constants,  namely,  and  K  coincide  with  the 

bounds,  whereas,  the  fifth,  which  is  the  transverse  shear  modulus,  Gj  ,  lies  in  between. 

b.  Stress  Field  Around  a  Fibrous  Inclusion 

Edwards  [49]  has  obtained  exact  closed  solutions  for  the  distribution  of  stress  around  a 
spheroidal  inclusion  in  an  elastic  body,  which  is  in  a  uniform  state  of  stress  infinitely  far  from  the 
inclusion.  For  the  special  case  of  uniaxial  tension  applied  perpendicular  to  the  major  axis  of  the 

inclusion,  as  shown  in  Fig.  28  below,  the  normal  stress  component,  ,  has  been  evaluated  at  the 
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two  equatorial  points,  B  and  C ,  as  a  function  of  the  parameter  H. 


By  definition. 


where 


e.g. 


H  =  -  1 

G 


G'  =  shear  modulus  of  the  inclusion 
G  =  shear  modulus  of  the  matrix 

H=  -1 ,  a  cavity 

H  =  0  ,  single  homogeneous  medium 
H  ->  oo ,  rigid  inclusion 


Table  5.2:  Comparison  of  stress  components  at  the  equator 
(Edward's  solution  and  present  analysis) 


0,atB 

(calculated) 

Gjj  at  B 

(exact  result) 

o,atC 

(calculated) 

CT,atC 

(exact  result) 

-  1.0 

3.15E-06 

- 

2.9999 

3.0 

-0.5 

0.7573 

0.743 

1.5058 

1.5 

0. 

1.0 

1.0 

1.0 

1.0 

10.0 

1.4158 

1.4 

0.1034 

0.089 

OO 

1.4770 

1.471 

-.0317 

-.035 

c.  Stress  Field  Around  a  Coated.  Continuous  Fiber 

The  stress  field  in  a  coated  continuous  fiber  composite  subjected  to  thermo-mechanical 
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loading  has  been  calculated  by  Mikata  and  Taya  [35]  by  use  of  four  concentric  circular  cylinder 
model.  The  target  material  was  Ni-  or  SiC-  coated  graphite  fiber/6061  aluminium  composite. 
From  the  material  properties  used,  it  is  noted  that  the  matrix  and  the  coating  are  isotropic,  but 
graphite  fiber  is  transversely  isotropic.  The  geometrical  model  used  in  their  analysis  is  illustrated  in 
Figure  29. 


d.  Uniaxial  Tension 

For  a  uniaxial  applied  stress,  ,  the  maximum  stresses  which  occur  at  point  A,  within 
the  coating,  have  been  evaluated  by  changing  two  parameters;  Vf  ,  the  volume  fraction  of  fiber,  and 

c/d,  which  is  the  ratio  of  coating  thickness  to  fiber  diameter.  Their  results  along  with  our  theoretical 
predictions  are  listed  in  table  5.3. 


Table  5.3:  Comparison  of  G/Cq  point  A  within  the  coating 
(Mikata  and  Taya's  results  and  present  analysis) 


Vf  c/d 

Ni  coating 

SiC  coating 

M  &T 

Our's 

M  &T 

Our's 

0.25 

1/70 

1.880 

1.872 

4.18 

4.141 

5/70 

1.750 

1.739 

3.41 

3.403 

0.36 

3/70 

1.550 

1.537 

3.14 

3.122 

5/70 

1.480 

1.468 

2.80 

2.784 

0.49 

3/70 

1.320 

1.308 

2.64 

2.610 

5/70 

1.260 

1.241 

2.31 

2.300 
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e.  Uniform temperaturechangef  AT >Q) 

As  a  special  case  of  the  fonnuladon,  Mikata  and  Taya  obtained  an  explicit  closed  form 
solution  for  a  single  inhomogeneity  problem  (i.e.  a  continuous  fiber  embedded  in  an  infinite  matrix). 
Using  the  material  properties  of  graphite  T300  fiber  and  A1  6061  matrix,  the  following  stress  field 

was  obtained  for  the  composite  subjected  to  a  uniform  temperature  change,  A  T. 


Table  5.4  :  Stress  comparison  for  a  uniform  temperature  change 
(Mikata  and  Taya's  results  and  present  analysis) 


Normalised 

Fiber 

Matrix  (at  interface) 

a/AT 

Oq/AT 

ojat 

a/AT 

GqIA 

a^AT 

M  &  T 

8.855 

8.855 

5.703 

8.855 

-8.855 

5.084 

E  +  4 

E  +  4 

E  +  6 

E  +  4 

E  +  4 

E- 10 

Present 

8.855 

8.855 

5.703 

8.855 

-8.855 

-1.425 

E  +  4 

E  +  4 

E  +  6 

E  +  4 

E  +  4 

E-02 

5.7  SUMMARY 

Generally  speaking,  very  good  agreement  has  been  obtained  between  our  results  and  the 

analysis  of  other  researchers  in  the  field.  What  is  specially  to  be  noted  is  the  fact  that  a  model  has 

been  developed  to  predict  the  effective  thermoelastic  properties  of  multidirectional  coated  fiber 

composites  and  the  solution  to  the  stress  distribution  has  been  obtained  under  a  uniform 
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three-dimensional  mechanical  and/or  non-mechanical  loading.  It  is  therefore  very  encouraging  to 
observe  that  the  analysis  does  reduce  down  to  the  available  exact  solutions  for  some  of  the  simpler 
cases.  This  has  therefore  given  us  reasonable  confidence  and  we  can  proceed  to  use  this  program  for 
the  analysis  and  parametric  study  of  fiber  reinforced  con^site  materials. 

5.8  THE  NDSANDS  PROGRAM 

Micromechanical  considerations  in  composite  materials  may  require  the  use  of  a  practical 
tool  that  can  handle  different  constituent  materials,  arbitrary  fiber  orientations  and  multi-axial  loading 
conditions.  To  address  these  requirements,  the  computer  code  called  NDSANDS  (H  Directional 
Stiffness  A  N  2  Strength)  has  been  developed.  A  flow  chart  of  the  program  is  given  in  Figure  30. 
As  seen  from  the  diagram,  the  program  can  be  used  either  to  analyze  a  composite  or  to  conduct 
parametric  study.  By  parametric  study  what  is  meant  is  that  the  user  can  change  either  a  material 
property  or  the  geometry  of  the  composite,  CMie  single  variable  at  a  time  while  the  remainder 
are  kept  constant,  and  thereby  examine  the  change  in  effective  properties  and  stress  distribution  as  a 
result  of  different  input  values  of  the  parameter  selected.  When  changing  the  material  property,  we 
must  insure  that  both  the  stiffness  and  compliance  matrices  remain  positive  definite  at  all  times  [50]. 

Description  of  the  Computer  Program 

The  computer  code  NDSANDS  was  modified  sc*  that  it  could  be  implemented  on  a 
IBM  PC  and  generalised  so  that  it  became  more  efficient,  user  friendly  and  easily  menu  driven. 
What  follows  is  a  brief  description  of  some  of  the  procedures  that  were  incorporated  in  the  original 
program. 

i)  Data  files  FIBRE,  SHEATH,  MATRIX  and  PARTICLE  which  enlisted  the  complete 
set  of  property  characterization  for  each  one  of  the  constituents,  were  created.  This  data  is  intended 
to  be  of  value  primarily  in  the  preliminaryt  selection  of  materials  to  fulfill  the  requirements  of 
anticipated  aerospace  applications.  A  comprehensive  list  of  the  constituent  properties  is  given  in  the 
next  section. 
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ii)  By  a  suitable  choice  of  constituent  material  from  each  one  of  the  material  libraries,  and 
by  specifying  the  geometric  variables,  such  as  radii,  orientation  angles  etc.,  the  user  can  create 
his/her  own  composite.  At  present,  one  can  choose  two  different  types  of  fibres  and  two  different 
types  of  sheath  materials  and  specify  a  maximum  of  ten  different  orientations  of  the  flbo'. 

iii)  A  maximum  of  five  different  composites  can  be  analyzed  simultaneously.  Their 
effective  thermoelastic  properties  are  calculated  and  presented  in  a  tabular  form  so  that  direct 
comparisons  can  be  made.  The  stress  distribution  under  an  externally  applied  general 
three-dimensional  mechanical  loading  and/or  non-mechanical  loading  (specify  AT,  the  temperature 
difference,  and  thereby  calculate  the  thermal  load  or  specify  AC,  the  concentration  difference,  and 
calculate  the  hygrothermal  load)  could  then  be  determined. 

iv)  The  stress  analysis  can  be  carried  out  either  in  cylindrical  or  cartesian  coordinates 
depending  upon  the  user  input.  The  computed  stress  components  can  be  listed  in  a  table  or  can  be 
plotted.  It  is  upto  the  user  to  define  the  grid  points  for  listing  the  stress  components  or  give  the  plot 
specifications  if  he/she  is  interested  in  seeing  the  graphs.  At  present,  one  can  plot  either  at  a  constant 
radius  (with  variable  angle)  or  at  a  constant  angle  (with  variable  radius).  The  user  can  change  the 
number  of  grid  points  for  listing  or  the  plot  specifications  at  any  time  when  needed.  Using  the 
optimization  routine,  COMPLEX  (FUNK),  the  location  where  each  one  of  the  stress  component 
(cylindrical  or  cartesian)  exhibits  a  maximum  or  a  minimum  can  be  further  determined. 

v)  The  maximum  stress  criteria  has  been  incorporated  into  the  irxxlel  to  predict  the  initial 
or  first  failure  of  the  composites.  This  criteria  is  applied  by  calculating  the  strength/stress  ratio  for 
each  stress  component.  The  lowest  strength  ratio  controls  the  failure  of  the  composite.  Similar  to 
the  stress  components,  the  user  has  the  option  to  list,  optimize  or  plot  any  or  all  of  the  strength 
ratios. 


vi)  Once  the  effective  thermoelastic  properties  of  the  composite  have  been  calculated,  the 
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user  has  the  option  of  changing  mechanical,  thermal  or  hygrothermal  loads  without  having  to 
recalculate  the  effective  properties.  The  original  program  could  handle  only  one  loading  condition 
applied  to  a  single  composite,  whereas,  now  with  the  present  routine,  five  different  composites  can 
be  analysed  simultaneously  under  different  loading  conditions. 

vii)  Finally,  the  computer  code  NDSANDS  has  been  set  up  to  conduct  parametric  study. 
A  composite,  as  already  pointed  out,  is  created  by  choosing  constituent  materials  and  by  specifying 
the  geometric  variables.  In  this  part  of  the  study,  the  user  can  change  either  the  material  property  or 
the  geometry  of  the  composite,  and  analyse  the  change  in  effective  properties  and  stress  distribution 
by  inputting  five  different  values  for  the  parameter  selected.  We  have  also  modified  the  code  to  plot 
the  effective  properties  of  the  composite  selected  for  analysis  as  a  function  of  the  parameter  under 
consideration  so  that  the  influence  of  the  variable  could  be  easily  discerned. 

5.9  NUMERICAL  RESULTS  AND  DISCUSSION 

The  second  phase  of  the  project  was  associated  with  the  actual  analysis  of  some  of  the 
advanced  composite  materials  available.  We  were  specially  interested  in  analysing  composites 
intended  for  high  temperature  applications. 

In  the  general  theory  developed  in  first  section  of  this  report,  the  fiber  and  the  coating  are 
transversely  isotropic  characterized  by  five  independent  elastic  constants,  whereas,  the  matrix  is 
isotropic  in  nature,  and  needs  only  two  constants  for  material  specification.  Besides  the  elastic 
properties,  we  are  also  interested  in  knowing  the  hygrothermal  properties  such  as  coefficient  of 
thermal  expansion,  moisture  expansion,  thermal  conductivity,  etc.,  of  the  constituents  to  predict  the 
effective  thermo-elastic  properties  of  the  composite.  Finally,  in  order  to  calculate  the  strength  ratios 
and  thereby  predict  the  initial  or  first  failure,  strength  characterization  in  the  longitudinal  and 
transverse  directions,  both  ,  under  tension  and  compression,  as  well  as  the  behaviour  in  shear,  is 
required  to  be  known,  for  all  the  concerned  phases. 
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Thus,  for  a  transversely  isotropic  constituent,  we  need  a  K>tal  of  sixteen  properties  to 
carry  out  the  study.  An  extensive  literature  survey  was  dierefore  undertaken  to  obtain  a  conq^lete  set 
of  data  for  each  phase.  One  of  the  major  problems  which  we  encountered  while  canying  out  this 
survey  was  how  to  compile  the  elastic,  thermal  and  strength  properties  of  different  constituents.  In 
most  of  the  reported  literature,  studies  have  been  undertaken  on  only  one  aspect  of  the  problem,  say 
the  elastic  property,  and  the  remaining  two  have  eidier  beat  ignored  or  not  mentioned  at  all. 

What  follows  now  is  the  data  that  we  were  able  to  get  together  for  different  phases  after 
an  exhaustive  search  [51  -  77].  It  should  be  pointed  out  that  it  is  still  not  complete.  This  results  from 
a  lack  of  extensive  testing  of  a  given  material  and  also  from  the  sensitivity  of  the  properties  of 
different  materials  to  fabricaticm  and  test  techniques,  so  that  comparism  of  data  from  various  sources 
was  difficult  This  data  is  intended  to  be  of  value  primarily  in  the  preliminary  selection  of  materials. 
It  has  been  compiled  to  assist  the  designer  in  selecting  and  specifying  one  or  more  material  that 
appears  to  fulfill  the  requirements  of  anticipated  aerospace  applicadtxis. 
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5.5  PROPERTY  TABLES 


FIBER  TABLE  -  1 


FIBER  NAME 

PROPERTIES 

:  T300(Gr) 

Modmw  II 

AS4(Carbon) 

Boron 

ET(Psi) 

:  0.3I90E407 

0.2170E+07 

0.2030E407 

0.6547E+08 

EA(Psi) 

0.3277E+08 

0.3370E+08 

0.3408E+08 

0.5405E+08 

NUT 

0.4200E+(K) 

0.4900E+00 

0.2500E400 

0.1210E-fO0 

NUA 

:  0.3000E+00 

0.2790E400 

0.2000E400 

0.1300E400 

GA(Psi) 

:  0.7003E406 

0.3480E407 

0.4008E407 

0.2742E408 

ALFAT(I/OF) 

:  0.1500E-04 

0.5000E-05 

O.lOOOE-04 

0.2700E-05 

ALFAA(l/(^0 

:  -0.8340E-06 

-0.2000E-06 

-0.2000E-06 

0.2700E-05 

BETAT(I/%M) 

:  O.OOOOE+00 

O.OOOOE+00 

O.OOOOE-fOO 

O.OOOOE400 

BETAA(I/%M) 

:  0.0000E-f00 

O.OOOOE400 

O.OOOOE+00 

O.OOOOE+00 

MUT(Btu/hr-ft-OF) 

0.483  lE+01 

0.3750E+01 

MUA(Btu/hr-ft-OF) 

0.483 1E402 

0.1386E+O2 

LULT(Psi) 

0.4408E+06 

0.4250E+06 

0.5207E+06 

0.460(p406 

LULC(Psi) 

0.3335E-f06 

0.3335E+06 

0.6000E+06 

TULT(Psi) 

0.4500E+05 

0.5076E+05 

0.3300E+06 

TULC(Psi) 

0.3000E-f05 

LS(Psi) 

0.6000E+04 

0.1200E+05 
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FIBER  TABLE -2 


FIBER  NAME 

PROPERTIES 

E-Giass 

AS  (Graphite) 

S2-Glass 

Beryllium 

ET(Psi) 

;  0.1053E+08 

0.2000E407 

0.1250E408 

0.4200E408 

EA(Psi) 

:  0.10S3E-f08 

0.3200E+08 

0.12S0E4O8 

0.4200E408 

NUT 

:  0.2200E+00 

0.2500E400 

0.220OE-tO0 

0.3000E-01 

NUA 

:  0.2200E+00 

0.2000E+00 

0.2200E+00 

0.3000E-OI 

GA(Psi) 

:  0.4316E-f07 

0.5000E407 

0.5100E-f07 

0.2040E408 

ALFAT(1/0F) 

:  0.2780E-05 

0.5000E-05 

0.2780E-05 

0.6450E-05 

ALFAA(1/0F) 

:  0.2780E-05 

-0:2000E-06 

0.2780E-05 

0.6450E-05 

BETAT(1/%M) 

:  O.OOOOE+00 

0.0(X)0E-fO0 

O.OOOOE-»4)0 

O.OOOOE-tOO 

BETAA(1/%M) 

:  O.OOOOE400 

O.OOOOE+00 

O.OOOOE-i4)0 

O.OOOOE+00 

MUT(Btu/hr-ft-OF) 

:  0.6250E400 

0.483 1E401 

0.4830E400 

0.8455E402 

MUA(Btu/hr-ft-OF) 

:  0.6250E+00 

0.483  lE+02 

0.4830E+00 

0.8455E+02 

LULT(Psi) 

0.4979E+06 

0.4500E+06 

0.7000E+06 

0.1300E-K)6 

LULC(Psi) 

;  0.7250E405 

0.3335E406 

0.7250E-fO5 

0.1400E406 

TULT(Psi) 

0.4979E406 

0.5000E405 

0.7000E+06 

0.1300E406 

TULC(Psi) 

LS(Psi) 

0.7250E+05 

0.2500E406 

0.2500E405 

0.7250E+05 

0.3500E-f06 

0.1400E+06 
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nBER  TABLE -3 


RBERNAME 

PROPERTIES 

Nicalon 

AVCO 

HM(Cart)on) 

ETCPsi) 

:  0.2900E+08 

0.6000E408 

0.1500E+07 

EA(Psi) 

:  0.2900E+08 

0.6000E408 

0.5200E-f08 

NUT 

:  0.2987E+00 

0.2500E+00 

0.3600E+00 

NUA 

:  0.2987E+00 

0.2500E+00 

0.2600E400 

GA(Psi) 

:  0.1117E+08 

0.2400E+08 

0.2090E+07 

ALFAT(1/0F) 

:  0.1778E-05 

0.2330E-05 

0.4390E-05 

ALFAA(1/0F) 

0.1778E-05 

0.2330E-05 

-0.4230E-06 

BETAT(1/%M) 

:  O.OOOOE+00 

O.OOOOE+00 

O.OOOOE-fOO 

BETAA(1/%M) 

:  O.OOOOE+00 

0.0000E4O0 

O.OOOOE-fOO 

MUT(Btu/hr-ft-OF) 

;  0.6710E401 

0.2367E+02 

MUA(Btu/hr-ft-OF) 

0.6710E+01 

0.2367E402 

LULT(Psi) 

:  0.4000E+06 

0.5000E+06 

0.2700E406 

LULC(Psi) 

:  0.2560E+06 

0.2559E-f06 

0.2900E-f06 

TULT(Psi) 

0.4000E-^06 

0.5000E406 

TULC(Psi) 

0.2560E+06 

0.2559E+06 

LS(Psi) 
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SHEATH  TABLE 


SHEATH  NAME 

PROPERTIES 

Ni 

SiC 

Polyimide 

Carbon 

ET(Psi) 

0.3001E408 

0.7003E+08 

0.4522E+06 

0.1320E+07 

EA(Psi) 

:  0.3001E+08 

0.7(X)3E408 

0.4S22E-f06 

0.1320E+07 

NUT 

:  0.3100E-fO0 

0.1900E400 

0.3300E-fO0 

O.llOOE+00 

NUA 

:  0.3100E-fOO 

0.1900E400 

0.330OE-fO0 

O.llOOE-fOO 

GA(Psi) 

:  0.1145E+08 

0.2943E408 

0.1700E-f06 

0.5946E-H)6 

ALFAT(1/0F) 

:  0.7390E-05 

0.241  lE-OS 

0.3000E-04 

0.1220E-05 

ALFAA(1/0F) 

:  0.7390E-05 

0.241  lE-05 

0.3000E-04 

0.1220E-05 

BETAT(1/%M) 

:  O.OOOOE-fOO 

0.0000E4O0 

BETAA(1/%M) 

O.OOOOE-fOO 

0.0000E4O0 

MUT(Btu/hr-ft-OF) 

0.1700E+01 

0.1200E403 

MUA(Btu/hr-ft-OF) 

0.1700E+01 

0.1200E+03 

LULT(Psi) 

:  0.4596E+05 

0.8004E40S 

0.5600E-t04 

LULC(Psi) 

0.7600E-^04 

TULT(Psi) 

:  0.4596E+05 

0.8004E+05 

0.5600E-f04 

TULC(Psi) 

0.7600E404 

LS(Psi) 

0.2920E404 
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MATRIX  TABLE 


MATRIX  NAME 

PROPERTIES 

EPON  828 

H3501-6 

X7002-T6(A1) 

Graphite 

E(Psi) 

:  0.4200E406 

0.6200E406 

0.1037E+08 

0.1410E407 

G(Psi) 

:  0.1560E+06 

0.2300E+06 

0.3900E+07 

0.6350E+06 

ALFA(I/0F) 

:  0.2667E-04 

0.2272E-04 

0.1334E-04 

0.3340E-05 

BETA(I/%M) 

:  0.3200E-02 

0.3200E-02 

MU(Btu/hr-ft-OF) 

0.1420E400 

0.1420E+00 

0.1280E+03 

0.1200E+03 

ULT(Psi) 

:  0.1025E+05 

0.1200E+05 

0.5400E-f05 

0.3000E404 

ULC(Psi) 

:  0.2610E+05 

0.300(ffi+05 

0.6750E+05 

0.7600E+04 

S(Psi) 

;  0.5125E+04 

0.6000E+04 

0.3240E+05 

0.3000E+04 

MATRIX  NAME 

PROPERTIES 

:  A  277-I5(P) 

Al  6061 

SiC 

E-Glass 

E(Psi) 

0.7330E+06 

0.9947E+07 

0.5600E+08 

0.1050E+08 

G(Psi) 

:  0.3330E+06 

0.3698E+07 

0.244  lE+08 

0.4303E+07 

ALFA(1/0F) 

0.1305E-04 

0.2500E-05 

0.1583E-04 

BETA  (1/%M) 

MU(Btu/hr-ft-OF) 

:  0.1200E+03 

0.2367E+02 

0.2050E400 

ULT(Psi) 

0.3330E+04 

0.4495E+05 

0.2250E+05 

0.1150E+05 

ULC(Psi) 

:  0.6670E+04 

0.8200E+05 

0.3400E+05 

S(Psi) 

0.6010E404 

0.3000E+05 

0.5750E+04 
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MATRIX  NAME 


:  ATJS(Carbon) 


Ccr-VIT 


Zinc  Bromide 


A1  Oxide 


PROPERTIES _ 

E(Psi) 

;  0.1320E407 

0.1340E-f08 

0.3632E-f08 

0.4500E+08 

G(Psi) 

:  0.5946E+06 

0.5317E467 

0.1600E408 

0.1800E+08 

ALFA(1/0F) 

:  0.1220E-05 

0.3300E-05 

BETA  (1/%M) 

MU(Btu/hr-ft-OF) 

:  0.1200E403 

0.%70E+00 

0.1500E-f02 

0.1500E402 

ULT(Psi) 

:  0.S600E-fO4 

0.1900E+06 

0.2250E+05 

0.4000E+05 

ULC(Psi) 

0.7600E-f04 

0.2000E-f06 

0.3500E+06 

S(Psi) 

:  0.2920E+04 

0.6500E405 

0.4000E+05 

MATRIX  NAME 

PROPERTIES 

Titanium 

Mg  Oxide 

LAS  ffl 

PMAS  III 

E(Psi) 

O.I500E+08 

0.2832E408 

0.1276E-f08 

0.1537E+08 

G(Psi) 

:  0.5770E407 

0.1200E408 

0.5221E+07 

0.6236E+07 

ALFA(1/0F) 

0.5000E-05 

0.8334E-06 

0.1500E-05 

BETA  (1/%M) 

MU(Btu/hr-ft-OF) 

:  0.9660E401 

0.2500E+02 

ULT(Psi) 

:  0.I350E+O6 

0.1400E+05 

ULC(Psi) 

0.2000E406 

S(Psi) 

: 

0.2400E405 
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MATRIX  NAME 

PROPERTIES 

:  1723  Glass 

MAS-LZN 

7052  Glass 

7050  Glass 

E(Psi) 

0.1276E408 

0.1500E408 

0.8200E+07 

0.8700E+07 

G(Psi) 

:  0.5221E+07 

0.6250E+07 

0.336  lE+07 

0.3566E+07 

ALFA(I/0F) 

BETA  (1/%M) 

MU(Btii/hr-ft-OF) 

ULT(Psi) 

ULC(Psi) 

S(Psi) 

:  0.2889E-05 

0.4300E-05 

0.5000E+04 

0.1000E-f06 

0.2872E-05 

0.2%7E-05 

MATRIX  NAME 

PROPERTIES 

:  7070  Glass 

7740  Glass 

1720  Glass 

9741  Glass 

E(Psi) 

:  0.7400E+07 

0.9100E+07 

0.1270E+08 

0.7200E+07 

G(Psi) 

:  0.3033E+07 

0.3792E+07 

0.5121E+07 

0.2927E+07 

ALFA(1/0F) 

0.2167E-05 

0.1944E-05 

0.2789E-05 

0.2722E-05 

BETA  (1/%M) 
MU(BUi/hr-ft-OF) 
ULT(Psi) 
ULC(Psi) 

S(Psi) 
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MATRIX  NAME 


CAS  I 


pboferhes 


E(Psi) 

:  0.1276E408 

G(Psi) 

:  0.5221E407 

ALFA(1/0F) 

:  0.2777E-05 

BETA  (1/%M) 

MU(Btu/hr-ft-OF) 

ULT(Psi) 

ULC(Psi) 

S(Psi) 

PARTICLE  TABLE 


PARTICLE  NAME 

WC 

Glass 

Quartz  sand 

E(Psi) 

:  0.1020E+09 

0.1020E408 

0.1067E-»08 

G(Psi) 

;  0.4180E408 

0.4215E+07 

0.4269E-M)7 

ALFA(1/0F) 

BETA(1/%M) 

MU(Btu/hr-ft-OF) 

ULT(Psi) 

ULC(Psi) 

S(Psi) 
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5.10  ILLUSTRATIVE  PROBLEM 


We  now  consider  a  three-dimensional  fibrous  composite  by  arranging  six  fibers  parallel 
to  the  six  lines  joining  the  q)posite  vertices  of  a  regular  icosahedron.  These  six  axes  can  be  oriented 

with  reference  to  an  orthogonal  cartesian  coordinate  system  x^  X2  X3  as  follows:  one  pair  in  the  X1X2 

plane  making  angles  of  6  with  the  x^-axis,  one  pair  in  the  X2  X3  plane  making  angles  of  6  with 

the  X2-axis,  and  one  pair  in  the  X3X1  plane  making  angles  of  6  with  the  X3-axis,  where  6 
=  tan-H2sin  180)  =  3l0  43'. 


As  shown  by  Rosen  and  Shu  [77],  this  type  of  arrangement  gives  rise  to  local  isotropy. 

The  isotropic  relation  G  =  E  /  [2.0  (1.0  +  v )]  can  be  used  as  an  independent  check  of  the  model.  In 
general,  this  relation  is  not  satisfied  exactly  in  the  present  analysis,  however  the  error  is  very  small. 
The  results  for  the  Nicalon  fiber  and  BMAS  matrix  system  arc  shown  in  table  5.6.  Also  shown  in  the 
table  is  the  effect  of  different  coating  thicknesses  and  coating  materials  on  the  effective  thermoelastic 
moduli.  The  fiber  volume  fraction  was  set  at  30%  but  the  ratio  of  coating  thickness  to  the  cylinder 

outer  radius,  defined  as  ( r2  -  r^ )  /r3 ,  was  treated  as  an  independant  variable. 

As  seen  from  table  5.6  below,  both  carbon  and  polyimide  coatings,  because  of  their 
lower  elastic  moduli  compared  to  the  BMAS  matrix,  reduce  the  effective  elastic  properties,  whereas, 
nick:*.!,  with  a  higher  moduli,  adds  to  the  reinforcement.  The  thermal  expansion  coefficient  of  the 
composite,  on  the  other  hand,  is  influenced  more  by  nickel  and  polyimide  coatings  because  of  the 
larger  degree  of  mismatch  between  the  constituents. 

The  microstress  distribution  within  the  constituents,  for  a  multidirectional  fiber 

composite,  in  general,  depends  both  on  the  type  of  loading  and  the  fiber  orientation.  As  an 

approximation,  the  curing  or  residual  stresses  can  be  estimated  by  subjecting  the  composite  to  a 

89 


uniform  temperature  change.  For  the  specific  three-dimensional  composite  under  consideration, 
the  stress  distribution  remains  identical  for  the  six  fiber  orientations  for  a  unit  degree  cool  down. 


In  this  problem,  the  only  non-zero  stress  components  predicted  by  the  present  riKxiel  are  ,  Cq  and 

Oy  .  The  stress  concentration  is  a  function  of  both  Young's  modulus  and  the  thermal  expansion 
coefficient  of  the  coating,  besides  its  thickness.  Hie  effect  of  different  coating  materials  and/or 
thicknesses  on  the  stress  concentrations  is  quite  dramatic.  These  trends  are  illustrated  in  tables 

5.7  (a),  5.7  (b)  and  5.7  (c)  for  the  components  c, ,  C0  and  ,  respectively. 

The  radial  stress  component  at  the  interface  can  be  considered  as  a  failure  criteria  for 

debonding,  e.g.,  a  negative  value  of  promotes  contact  between  the  constituents,  whereas,  a 

positive  value  suggests  possible  separation  and  initiation  of  debonding  at  the  boundary.  It  is  seen  that 
a  'thick'  coating  of  a  'soft'  material  with  a  'low'  coefficient  of  thermal  expansion  helps  in  reducing 
the  stress  concentration  factor  at  the  boundary.  Within  the  coating,  the  algebraic  maximum  hoop 
stress  occurs  at  the  fiber-coating  interface,  whereas,  in  the  noatrix,  the  maximum  occurs  at  the 
coating-matrix  interface.  Even  though  the  polyimide  coating  has  a  low  elastic  modulus  (3.1  GPa) , 
but  because  of  its  high  coefficient  of  thermal  expansion  (  54  x  10'^/ ),  the  magnitude  of  the 
constituents  stress  components  for  the  Nicalon/Polyimide/BMAS  composite  are  still  large  in 
magnitude  for  a  unit  degree  cool  down. 


To  conclude,  it  is  apparent  that  generally  a  reduction  in  the  stress  concentration  can  be 
made  at  the  expense  of  the  elastic  moduli  of  the  composite.  Further,  by  a  proper  choice  of  coating 
thickness,  modulus  and  coefficient  of  thermal  expansion,  the  stress  component  of  interest,  which  is 
instrumental  in  causing  a  specific  mode  of  failure,  can  be  controlled. 
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TABLE  5.6:  Three>(iimensional  '  isotropic  '  composite,  effective  moduli 


Coatine  thickness 

Cylinder  outer  radius 

Composite  system 

E 

(GPa) 

V 

G 

(GPa) 

a 

(10-6/OC) 

0 

Nicalon/BMAS 

128.16 

0.250 

51.25 

2.89 

Nicalon/Nickel/ 

BMAS 

129.11 

0.251 

51.60 

3.04 

0.01 

Nicalon/Caiton/ 

BMAS 

118.06 

0.243 

47.50 

2.88 

Nicalon/Polyimide/ 

BMAS 

106.01 

0.249 

42.45 

3.24 

Nicalon/Nickel/ 

BMAS 

138.72 

0.259 

55.10 

4.45 

0.10 

Nicalon/Carbon/ 

BMAS 

80.20 

0.217 

32.94 

2.84 

Nicalon/Polyimide/ 

BMAS 

65.93 

0.234 

26.72 

4.76 
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TABLE  5.7  (a)  :  Stress  component  in  the  fiber-coating  and  coating-matrix 
interface  for  AT  =  -1®C 


coating  thickness 

cylinder  outer  radius 

Conqwsite  system 

(Or)f^ 

(KPa) 

(KPa) 

0 

Nicalon/BMAS 

26.4 

26.4 

Nicalon/Nickel/ 

BMAS 

-8.98 

44.9 

0.01 

Nicalon/Carbon/ 

BMAS 

22.8 

22.3 

Nicalon/Polyimide/ 

BMAS 

104.0 

107.0 

Nicalon/Nickel/ 

BMAS 

-290.0 

140.0 

0.10 

Nicalon/Carbon/ 

BMAS 

8.54 

5.95 

Nicalon/Polyimide/ 

BMAS 

293.0 

306.0 
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TABLE  5.7  (b)  :  Stress  component  Cq  (algebraic  maximum)  in  the  fiber,  coating 


and  matrix  for  AT  =  -1®  C 


coadne  thickness 

cylinder  outer  radius 

Composite  system 

(O0)f 

(KPa) 

(Cf0)c 

(KPa) 

(<ye)m 

(KPa) 

0 

Nicalon/BMAS 

26.4 

- 

-55.6 

Nicalon/Nickel/ 

BMAS 

-8.98 

3020.0 

-97.5 

0.01 

Nicalon/Carbon/ 

BMAS 

22.8 

-7.46 

-50.7 

Nicalon/Polyimide/ 

BMAS 

104.0 

289.0 

-163.0 

Nicalon/Nickel/ 

BMAS 

-290.0 

2720.0 

-407.0 

0.10 

Nicalon/Caiton/ 

BMAS 

8.54 

-9.63 

-32.0 

Nicalon/Polyimide/ 

BMAS 

293.0 

383.0 

-494.0 
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TABLE  5.7  (c):  Stress  component  a,  in  the  fiber,  coating  and  matrix  for  AT  = 


coating _ thickness 

cylinder  outer  radius 

Composite  system 

(Oz)f 

(KPa) 

(Oz)c 

(KPa) 

(KPa) 

0 

Nicalon/BMAS 

77.9 

- 

-26.9 

Nicalon/Nlckel/ 

DMAS 

26.3 

3060.0 

-48.4 

0.01 

Nicalon/Carbon/ 

BMAS 

77.9 

-4.52 

-25.5 

Nicalon/Polyimide/ 

BMAS 

54.9 

288.0 

-69.7 

Nicalon/Nickel/ 

BMAS 

-424.0 

2590.0 

-248.0 

0.10 

Nicalon/Cartwn/ 

BMAS 

76.9 

-5.99 

-21.0 

Nicalon/Polyimide/ 

BMAS 

-136.0 

377.0 

-262.0 
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5.11  CONCLUDING  REMARKS 


In  summary,  we  have  developed  a  first  order  ideal  material  model  to  approximate  the 
elastic  response  of  a  composite  body  reinforced  by  coated  fibers  oriented  in  various  directions.  The 
coating  can  either  be  applied  intentionally  to  achieve  the  desirable  composite  pr(^)erties  or  it  can  occur 
due  to  the  processing  conditions  involved  in  the  composite  manufacture.  We  have  further 
demonstrated,  through  a  parametric  study,  how  an  applied  coating  modifies  the  micro  stress 
distribution  and  the  elastic  properties  of  a  three-dimensional  Nicalon/BMAS  system.  The  model  can 
be  used  to  provide  material  guidance  to  control  the  micromechanicul  failure  mode.s.  In  conjunction 
with  a  disciplined  experimental  program,  studies  such  as  those  conducted  here,  can  be  employed  to 
direct  further  improvements  in  the  model,  such  as  the  capability  to  handle  discontinuity  of  some  of 
the  traction  and  /or  displacement  components  at  the  boundary.  This  formulation  wiU  be  the  subject  of 
future  work  to  be  undertaken. 
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MAIN  PROGRAM 


INPUT  SUBROUTINE 


STIFFNESS  SUBROUTINI 


END  of  PROGRAM 


1  initial  program  data 
and  read  element  data 

2  CALL  MATPRO  SUBROUTINE 
:  input  material  property 

3  CALL  GAUSPT  SUBROUTINE 
input  GAUSS  points 

CALLGLOBAQ 

1  CALL  DMATPS- 

^  CALL  LOCAQ 


input  G/ 


_ 

t _ 

FRONTAL  SUBROUTINE 

2  CALL  SFR2  :  input  shape  function 

3  CALL  JACOB  :  calculating 

transformation  matrix 

1  if  body  forces  exist  then  CALL  BODYF 

2  if  transversed  loads  exist  then 

CALL  LOADPB  subroutine 

3  if  edge  loads  exist  then  CALL  EDGEF 

1  pre-frontal  data  and  boundary  code 

are  coded  first 

2  standard  frontal  procedure  is  coded 
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DEFINITION  OF  THE  NOTATION 


NBA.  NRB 

input  files  unit 

NFOA.NFOB  : 
NFOC,  NFOD 

output  files  unit 

MELEM 

maximum  element  number  per  layer 

MPOIN 

maximum  node  numbo*  per  layer 

MPLY 

total  number  of  plies 

MDIME 

maximum  number  of  dimension 

NDOFP 

maximum  degree  of  freedom  per  node  at  local  region. 

MGASP 

maximum  number  of  Gauss  Integral  points 

MEVAB 

maximum  dimension  of  stiffness  matrix  in  element  level 

MFRON 

maximum  frontal  size 

MTOTV 

total  maximum  variable  number 

MNODE 

maximum  number  of  nodes  per  element 

DATA  FROM  UNIT  =  NHA 


BLOCK  A 
VARIABLE  = 

XL  XY  NX  NY  NARIDX  NGRZDY  IDOL 


This  block  read  grid  data  and  number  of  plies  of  the  structure. 


XL 

YL 

NX 

NY 

NGRIDX 

NGRIDY 

IDGL 

KSTRAN 


X  coordinate 
Y  coordinate 

number  of  nodes  per  side  in  X  direction 
number  of  grids  per  side  in  Y  direction 
number  of  grid  in  X  direction 
number  of  grid  in  Y  direction 
total  number  of  plies  of  the  structure 


I 


KSTRAN 
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BLOCK  B 
VARIABLE  = 

JLIPA  INCT 

This  block  data  use  to  determine  the  shape  function  is  linear  or  parabolic  function. 

JLIPA  :  1  linear  shape  function 

2  quadratic  shape  function 

INCT  increment 

BLOCK  C 
VARIABLE  = 

NVnX  NCASE  NTYPE  NNODE  NDOFN  NMATS 

NPROP  NGAUS  NDIME  NSTRE 

This  block  read  some  control  data  including  material,  dimension,  ingeration,  boundary  and 
load  control  data. 


NVFDC  : 

total  number  of  boundary  points 

NCASE 

number  of  load  case 

NTYPE  : 

parameter  for  output 

NNODE  : 

number  of  nodes  per  element 

NDOFN 

total  degree  of  freedom  per  element 

NMATS 

parameter  for  different  material  region 

NPROP 

number  of  material  properties 

NGAUS 

number  of  gauss  point 

NDIME 

number  of  dimension 

NSTRE 

dimension  of  stress  vector 
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BLOCK  D 
VARIABLE  = 

NUMEL  MATNO  LNDDS 

This  block  read  material  ID  and  element  connectivity. 

NUMEL  element  number 

MATNO  region  ID 

LNODS  element  connectivity  vector 

BLOCK  E 
VARIABLE  = 

I  LPSN 

This  block  read  element  ID  and  element  connectivity  ftM*  linear  case 

I  element  ID 

LPSN  element  connectivity  vector 

BLOCK  F 
VARIABLE  = 

IPOIN  COORD 

This  block  read  nodes  coordinates 

IPOIN  node  ID 

COORD  node  coordinates  vector 


I 
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BLOCK  G 
VARIABLE  = 

NOFIX  IFPRE  PRESE  KODE 

This  read  general  boundary  condition  code 
NOFIX  node  number 

IFPRE  prescribed  code 

1 :  prescribed  displacement  or  load 

2 :  not  prescribed 

PRESC  prescribed  value 

KODE  index  of  node 

1  :  1  degree  of  freedom  per  that  node 

3  :  3  degree  of  freedom  per  that  node 
6  :  6  degree  of  freedom  per  that  node 

BLOCK  H 
VARIABLE  = 

NOEND  KJEND  ENDPS  KCOD 

This  block  read  boundary  code  for  last  plies  or  last  layer  of  symetry  structure 
NOEND  index  number 

KJEND  prescribed  code 

1  :  prescribed  displacement  or  load 

0 :  not  prescribed 

ENDPS  prescribed  value  of  displacement  or  load 

KCOD  index  of  node 

1  :  1  degree  of  freedom  of  that  node 

3  :  3  degree  of  freedom  of  that  node 

6  :  6  degree  of  freedom  of  that  node 
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DATA  FROM  UNIT  =  NFIB 

BLOCK  I 
VARIABLE  = 

NPLY  NGROPG  NGROPL  NINDEX 

This  blcKk  read  group  data 
NPLY  total  number  of  group 

NGROPG  :  total  number  of  global  group 
NGROPL  :  total  number  of  local  group 
NINDEX  : 

BLOCK  J 
VARIABLE  = 

NGS  NGE  NLSE 

This  block  read  control  data  for  global  and  local  region 
NGS  starting  layer  number  of  global  vector 

NGE  ending  lay  number  of  global  vector 

NLSE  starting  layer  number  of  local  vector 

BLOCK  K 
VARIABLE  = 

JPLY  IGLCOD  PROPS 

This  bolock  read  ply  geometry  and  properties  • 

JPLY  :  ply  ID 
IGLCOD  region  ID 

PROPS  material  properties  vector  including  angle,  thickness,  modulus  poisson  ration 

shear  modulus 

I 
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Figure  1.  Global  and  Local  Regions  of  Laminate 


log 


=  Ai/A 

A  =  area  of  triangle 


Figure  2,  Triangular  (Area)  Coordinates 


Layer  Prop. 
Elem.  Geom. 


Moduli  for 

Global  Layer  GLOBQ 


Start  Element 
Calculations 


Obtain  Stiffness  Matrix  Global  Layer  -  VOGLOB 


Element  Layer  Assembly  -  ELSTG 


Material  Compliance  for  Local  Layer  k 


Stiffness  Matrix  Calculations  for  Local  Layer  k 


Element  Layer  Assembly  -  ELSTL 


Consistent  Load  Vector  -  LOAD 


Assembly  of  Element  Stiffness  Matrix  into 
Structure  Stiffness  Matrix  -  ASSMB 


Apply  Boundary  Conditions  -  BCOND 


Solution  of  Simultaneous  Equations  -  LEQTIP 


Listing  of  results  -  OUTPUT 


Figure  4.  Flow  Chart  for  Stiffness  Finite  Element  Method 
with  Global  and  Local  Layers 
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GRID  1  QRIO  2 


GRID  4 

Figure  5.  Finite  Element  Grid  Pattern  for  Square  Plate 
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(a)  n=2  (b)  n=4  (c)  n=6  (d)  n=8  (e)  n=10  (0  n=12 

(  n  :  number  of  elements  ) 


Figure  6.  Finite  Element  Grid  Pattern  for  Rectangular  Plate 
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n:  number  of  e 


Figure  7.  Convergence  Study  for  Stiffness 
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Figure  8.  Laminate  Geometry  for  Free  Edge 
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Figure  9.  Convergence  Study  of  Transverse  Displacement 
for  Rectangular  Plate 
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Max.  transverse  displacement 


Figure  11.  Convergence  Study  of  Normal  Stress  at  Mid-surface 

for  TensUe  Load 
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0.02-1 


Figure  12.  Convergence  Study  of  Normal  Stress  at  0/90  Interface 

for  Tensile  Load 
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Convergence  Study  of  Shear  Stress  in  z-y  Plane 
at  0/90  Interface  for  Tensile  Load 
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Figure  14.  Convergence  Study  of  Shear  Stress  in  z-x  Plane 


Weighted  displacement 


Figure  16.  Convergence  Study  of  Weighted  Displacement 
across  Top  Surface  for  Tensile  Lx)ad 
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Figure  17.  Comparing  Normal  Stress  (Mid-surface)  of  Different 
Geometry  Parameter  for  Tensile  Lx)ad 
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Figure  18.  Convergence  Study  of  Normal  Stress  at  Mid-surface 

for  Stretch  Load 


Figure  20.  Convergence  Study  of  Shear  Stress  in  z-y  Plane 
at  0/90  Interface  for  Stretch  Load 
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Figure  21.  Convergence  Study  of  Shear  Stress  in  z-x  Plane 
at  0/90  Interface  for  Stretch  Load 
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Weighted  displacement 


Figure  24.  Comparing  Normal  Stress  (at  0/90  Interface)  of  Different 
Geometry  Parameter  for  Stretch  Load 


Figure  25.  Comparing  Shear  Stress  (at  0/90  Interface)  of  Different 
Geometiy  Parameter  for  Stretch  Load 
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Figure  26.  Theoretical  Micro-mechanics  Model 
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Figure  27.  Composite  and  Equivalent  Body  Characterization 
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Figure  29.  Four  Riase  Model 
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Figure  30.  Flow  Chart  for  Micro-mechanics  Calculations 
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